**A Regression Model of Stream Water Quality Based on Interactions between Landscape Composition and Riparian Bu**ff**er Width in Small Catchments**

**Teresa Cristina Tarlé Pissarra 1,6, Carlos Alberto Valera 1,2,6, Renata Cristina Araújo Costa 1,6, Hygor Evangelista Siqueira 1,6, Marcílio Vieira Martins Filho 1,6, Renato Farias do Valle Júnior 3,6, Luís Filipe Sanches Fernandes 4,6 and Fernando António Leal Pacheco 5,6,\***


Received: 22 June 2019; Accepted: 20 August 2019; Published: 23 August 2019

**Abstract:** Riparian vegetation represents a protective barrier between human activities installed in catchments and capable of generating and exporting large amounts of contaminants, and stream water that is expected to keep quality overtime. This study explored the combined effect of landscape composition and buffer strip width (L) on stream water quality. The landscape composition was assessed by the forest (F) to agriculture (A) ratio (F/A), and the water quality by an index (IWQ) expressed as a function of physico-chemical parameters. The combined effect (F/A × L) was quantified by a multiple regression model with an interaction term. The study was carried out in eight catchments of Uberaba River Basin Environmental Protection Area, located in the state of Minas Gerais, Brazil, and characterized by very different F/A and L values. The results related to improved water quality (larger IWQ values) with increasing values of F/A and L, which were not surprising given the abundant similar reports widespread in the scientific literature. But the effect of F/A × L on IWQ was enlightening. The interaction between F/A and L reduced the range of L values required to sustain IWQ at a fair level by some 40%, which is remarkable. The interaction was related to the spatial distribution of infiltration capacity within the studied catchments. The high F/A catchments should comprise a larger number of infiltration patches, allowing a dominance of subsurface flow widespread within the soil layer, a condition that improves the probability of soil water to cross and interact with a buffer strip before reaching the stream. Conversely, the low F/A catchments are prone to the generation of an overland flow network, because the absence of permanent vegetation substantially reduces the number of infiltration patches. The overland flow network channelizes runoff and conveys the surface water into specific confluence points within the stream, reducing or even hampering an interaction with a buffer strip. Notwithstanding the interaction, the calculated L ranges (45–175 m) are much larger than the maximum width imposed by the Brazilian Forest Code (30 m), a result that deserves reflection.

**Keywords:** water pollution; riparian buffer width; landscape composition; regression model; interaction term; Brazilian Forest Code

#### **1. Introduction**

Riparian buffers represent undeniable benefits to stream water quality in catchments affected by agricultural nonpoint source pollution [1–4]. The benefits occur because riparian buffer strips favor physical processes, such as infiltration, sediment deposition or adsorption, as well as biochemical mechanisms, such as nutrient uptake or denitrification [5–7]. The role of riparian vegetation in the retention of nutrients and sediments has been reviewed in various studies [8–12]. A key parameter of riparian buffers is the width. Several studies refer minimum thresholds for riparian buffer width [13], while in some countries, such as Brazil or the United States, this width has been legally imposed or recommended [14,15]. However, other studies consider fixed-width recommendations problematic because riparian ecological responses are highly variable [8], and hence, optimal buffer widths are expected to vary from site to site [13]. On the other hand, width is just one among other major factors influencing stream water quality. Other key factors are landscape composition and patterns.

Following an early work by Kuehne, dated from the 1960s [16], numerous studies investigated the impacts of landscape composition and patterns on stream and lake water quality [17], even reporting that landscape characteristics are critical to water quality [18,19]. The reports evolved from cases where the relationship between the composition of landscape and the variation in water quality indicators was explored [20], to cases where the focus was moved to the spatial arrangement of the landscape (patterns) [21,22]. In the earlier studies, good water quality was generally associated with undeveloped watersheds dominated by forest land use, while poor water quality was linked to human development activities, such as agriculture [23]. In the more recent studies, a variety of landscape metrics was used to explain the correlation between landscape patterns and stream water quality, including patch density, largest patch index, landscape shape index or contagion [24–27].

The assessment of stream water quality based on the correlation with buffer strip widths or on the relationship with landscape composition usually leads to distinct results, and sometimes the standpoints are antagonistic. Various studies have shown that landscape composition within the river basin are decisive to identify the impacts of human activities on water quality [28,29], while other studies stated that patterns at the riparian buffer zones are more powerful to explain those impacts [30,31]. Eventually, the observed discrepancies were related to the dataset structure. Riparian variables are expected to become better explanatory variables when the land use is fairly homogenous and/or one land use category is widely dominant, as occurred in the [30,31] studies. In other cases, the landscape composition is almost always the first explanatory variable. However, discrepancies can also be interpreted as a scale problem [32]. The combined roles of the whole basin and buffer strip scales have been discussed in recent studies [33–36]. It has been reported that water quality varies along the direction of flow, due to human activities and the changing size of the buffer zone, and the self-purification ability of the water is influenced by the landscape composition in the river basin [37].

Despite the abundance of scientific literature on the relationship between stream water quality, buffer strip widths and landscape composition, a specific issue has not been tackled so far: The potential interaction among buffer strip widths and landscape composition. The studies mentioned in the previous paragraph link stream water quality variations to changes in buffer strip width and(or) landscape composition, acting independently as main effects, but fail to address potential joint effects. However, interactions among main effects are widely discussed in the scientific literature about regression models [38–40], and can play a role in the context of water quality assessments and their correlation with multi-scale factors. For example, stream water quality could be improved with narrower buffer strips if an enhanced self-purification of runoff was accomplished within a rural catchment by a large proportion of forest areas relative to agricultural areas. The negative consequences of overlooking interaction effects have been discussed in some forums (https://statisticsbyjim.com/regression/interaction-effects/). When interaction effects are statistically

significant, the main effects cannot be interpreted without considering the interactions, under the sentence of severe misinterpretation of results and prognosis.

The general purpose of this study was to explore the relationship between water quality variation, landscape composition, buffer strip widths, and potential interactions between composition and widths. This general goal comprised the following specific objectives: (1) Investigate potential interactions between landscape composition and buffer strip widths using a linear regression model with and without an interaction term; (2) determine the range of buffer strip widths that ensure a regular water quality, as a function of a landscape composition index (ratio native forest/agriculture), with and without considering the interaction effects; (3) interpret the results from a management standpoint. The research was carried out in the Uberaba River Basin (state of Minas Gerais, Brazil), namely within the Municipal Environmental Protection Area (EPA-MURB) located at the headwaters. The study area comprises eight sub-basins of EPA-MURB. Watercourses in these catchments may be affected by a diversity of non-point (diffuse) pollutants, including nitrogen and phosphorus from fertilizers or fine sediments from soil erosion, mostly derived from sugar cane plantations. In this study, water quality was assessed by an index that involves the measurement of dissolved oxygen, turbidity, total dissolved solids, which means parameters that can be interpreted as proxies to those pollutants. The index is called IWQ—Index for Water Quality and was proposed by the Environmental Company of São Paulo State—CETESB (https://cetesb.sp.gov.br) to be used in water quality assessments. The regression models were first applied to the IWQ and then to its formation variables, with the purpose to identify the most influencing ones. The studied watercourses are characterized by approximately 15, 30 and 50 m wide riparian forests. The reason for selecting these buffer widths relates to the rules imposed in the Brazilian Forest Code (Law No. 12651/12). There are two rules: The transition rule takes into account the size of land property calculated as fiscal modules and creates a distance from the stream margin that goes from a minimum of 5 m to a maximum of 20 m, considering the regular river bank; the permanent rule defines 30 m as unique distance. The rationale was, therefore, to span buffer widths that enclose these limits, considering the real buffer widths observed in the field.

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

#### *2.1. Study Area*

The study area encompasses the Municipal Environmental Protection Area of Uberaba River Basin (EPA-MURB), which covers 528.1 km2 of Triângulo Mineiro Region, State of Minas Gerais, Brazil (Figure 1). The EPA-MURB is located in the Brazilian Central Plateau, and the Northeast portion of Paraná River Basin, between the altitudes 710–1040 m and within the planimetric coordinates 188–220 km East and 7815–7840 km North expressed in the Universal Transverse Mercator system, zone 23K. The topography is characterized by undulated landscapes, sometimes incised by steep valleys.

Geology is dominated by two groups and associated formations (Figure 2a). The São Bento Group and associated Serra Geral Formation is composed of Lower Cretaceous grey to black basaltic lava flows (15–70 m thick), cropping out at lower altitudes. The Bauru Group, which overlays the São Bento Group along with an erosive contact, comprises the Uberaba Formation overlaid by the Marília Formation, being both composed of Cretaceous sandstones and conglomerates. The contact between the two sequences is abrupt, being marked by a silexite level and a conglomerate rich in quartz grains cemented with calcite [41]. The latosols dominate the landscape while the argisols occupy just small areas (Figure 2b) (https://www.embrapa.br/solos/sibcs). The latosols are characterized by clayey texture, whereas the argisols are characterized by sandy texture. Soil erosion may be intense because the land is prepared for seeding in the Autumn season, a period characterized by erosive rainfall events [42].

**Figure 1.** Location of the Environmental Protection Area of Uberaba River Basin in the Uberaba Municipality, State of Minas Gerais, and Brazil. Adapted from Valera et al. [15].

Climate is tropical (Aw in the Köppen's classification scheme) and the climatic domain is semi-humid with 4 to 5 dry months and relative humidity of 70–75%. The temperature ranges between 20 and 24 ◦C, on the annual average. The period from October to February observes the warmest temperatures that vary between 21 and 25 ◦C. The coldest month arrives in July, when temperatures drop to 16 to 18 ◦C. Based on a sixty-two years record, mean annual precipitation of 1584.2 mm is estimated for Uberaba municipality. The amount of rainfall varies considerably during the year, with average values ranging between 42.8 and 541 mm (www.inmet.gov.br/).

Agriculture and livestock production are important economic activities in the EPA-MURB. These areas form a mosaic where they alternate with spots of native vegetation (Cerrado), as illustrated in Figure 3. The landscape has changed significantly in the past half-century. In the 1960s, the Cerrado was the dominating land cover, reaching 40% of the EPA-MURB. In the following decades, expansion of managed pastures and (to a smaller extent) areas used for short-cycle agriculture (mostly corn and rice) caused a significant reduction in the share of native vegetation in the region. More recently, large areas have been converted to sugar cane plantations related to the production of energy from ethanol [43].

**Figure 2.** Geology (**a**) and soils (**b**) in the Environmental Protection Area of Uberaba River Basin. The soil map was adapted from Siqueira et al. [44].

**Figure 3.** Land cover in the Environmental Protection Area of Uberaba River Basin. Adapted from Siqueira et al. [44].

On 20 January 1999, the EPA-MURB was legally protected through the State Law No. 13183, for an area of 463 km2, being recognized as Sustainable Land Use Conservation Unit because of its important groundwater resources and remnants of Cerrado Biome. On 28 December 2005, the protection has been considered at Municipal level by the Law No. 9892 for an area of 528.1 km2. The management plan for the EPA-MURB was published by the Municipal Secretary of the Environment (http://www.uberaba.mg.gov.br/), which divided the area into five zones: (1) Urban consolidation zone; (2) tourism and/or leisure development zone; (3) agricultural area; (4) conservation area of natural resources and (5) recovery zone. Recently (2017), Complementary Law No. 561 has regulated the urban perimeter within the EPA-MURB.

#### *2.2. Studied Sub-Basins*

The studied sites comprised eight small to medium sub-basins selected within the EPA-MURB (Figure 4), with areas ranging from 136.3 hectares (Borá 1 sub-basin) and 3764 hectares (Lajeado sub-basin). The distribution of main land uses or covers is summarized in Table 1. In all cases the catchments were mostly used for agriculture, namely sugar cane plantations, as well as natural or managed (used for the grazing of domestic livestock) pastures. The use for agriculture ("A" column in Table 1) represents 56.7% (Borá 2 sub-basin) to 88.8% (Borá 2 sub-basin) of the catchment areas, and therefore, they can be considered basins with significant anthropogenic influence. Besides the use for agriculture, the eight catchments are substantially covered with native (Cerrado) and managed (eucalyptus) forests ("F" column), with proportions to the agriculture use ("F/A" column) ranging from 0.1 to 0.7. The other uses or covers ("Other" column) comprise the rural dwellings, water bodies (lakes and reservoirs), roads, orchards and areas used for rainfed or irrigated corps. The riparian buffers marginal to the watercourses ("L" column) were characterized by quite different widths: On average and approximately, 15 m in the Alegria, Lageado and Mangabeira 1 sub-basins; 30 m in the Borá 1, Mangabeira 2 and Uberaba sub-basins; and 50 m in the Borá 2 and Lanhoso sub-basins. The water samples were collected at the sub-basins' outlets.

**Figure 4.** Distribution of studied sub-basins within the Environmental Protection Area of Uberaba River Basin.

**Table 1.** Land use and cover within the studied sub-basins. Symbols: L—approximate average width of riparian buffer marginal to the water course, observed in satellite images; F—forest cover (native or managed); A—use for agriculture (sugar cane plantations; natural and managed pastures); Other—other uses or covers.


#### *2.3. Water Sampling and Physico-Chemical Analyses*

The stream water samples were collected in the streams, 60 cm far from the stream margin, in sectors that were adjacent to the riparian buffer. The sampling campaigns were conducted on a monthly basis, from January 2016 to January 2017 (13 months). The sampling took place from calendar days 15 and 20, every month. The year of 2016 was dry because the annual rainfall (1214.4 mm) was smaller than the local long-term average precipitation (1584.2 mm). Table 2 reviews the prevailing weather conditions in the sampling days, as well as during the three antecedent days. The sampling days were characterized by low rainfall (<5 mm), with the exception of February 2016 and January 2017 campaigns. In these two campaigns, rainfall during the sampling day has reached 5.9 and 10.9 mm, respectively. Average rainfall in the three antecedent days was also small (3.5–5.8 mm), with the few exceptions represented in Table 2 in boldface. The exceptional days were 13 November 2016, and 16 January 2017, with precipitation >25 mm. Taken altogether, the average analytical results should be related to long-term effects of landscape composition and buffer strip width on the quality of stream water, rather than short term effects associated with storm events.


**Table 2.** Weather conditions (rainfall) in the water sampling day and the three antecedent days (day-1 until day-3). Values larger than 10 mm day-1 are represented in boldface. Source: Valera et al. [15].

The measurement of water quality parameters in every campaign involved 10 repetitions, as recommended in the CONAMA's Resolution No. 357/2005. The parameters were measured using a Horiba U-50 Series multi-parameter probe, and comprised: T—water temperature (◦C), pH, ORP—Oxidation Reduction Potential (mV), Ec—Electrical conductivity (μS cm<sup>−</sup>1), Turbidity, measured in Nephelometric Turbidity Units (NTU), DO—Dissolved oxygen (mg L−1), PDO—Percentage of Dissolved Oxygen (%), and TDS—total dissolved solids (mg L<sup>−</sup>1). The analytical results are provided as Supplementary Material.

A subset of parameters was used to calculate the Index for Water Quality (IWQ) proposed by the Environmental Company of São Paulo State—CETESB (https://cetesb.sp.gov.br):

$$I\mathcal{W}\mathcal{Q} = \prod\_{i=1}^{n} q\_i^{w\_i} \tag{1}$$

where 0 ≤ *IWQ* ≤ 100, *qi* is the quality of *i* th parameter obtained from standardization of the measured values into a 0–100 range, *w*<sup>i</sup> is the weight of *i* th parameter, which varies in the 0 <sup>≤</sup> *w*<sup>i</sup> <sup>≤</sup> 1 interval as a function of its importance to the overall quality, and *n* is the total number of parameters. This *IWQ* evaluates water quality on the basis of nine parameters (*n* = 9), including water temperature, pH, dissolved oxygen, turbidity, total dissolved solids, biochemical oxygen demand, fecal coliforms, total nitrogen and total phosphorus. Data may be lacking for some parameters, but the index can still be calculated using the available data and adjusting the weights to different values. In the present study, the calculus of *IWQ* was based on five parameters (*n* = 5), namely the first five from the CETESB list, while the weights were set to—0.10 (water temperature); 0.21 (pH); 0.17 (turbidity); 0.2 (dissolved oxygen); 0.17 (total dissolved solids), as proposed in [45]. The transformation of measured parameters into *q* scores (Equation (1)) is based on standardization curves, which are portrayed in Figure 5 for the selected parameters. The average of each parameter in the studied catchments, as well as the corresponding *IWQ* values, are depicted in Table 3.

**Figure 5.** Standardization curves used to transform the water quality parameters into *q* scores (Equation (1)). Source: https://cetesb.sp.gov.br.


**Table 3.** Average water quality based on the five parameters used to calculate the IWQ (Index for Water Quality) index (Equation (1)). The full inventory of values, obtained within the monitoring period (January 2016–January 2017), is provided as Supplementary Material.

According to https://cetesb.sp.gov.br, the quality of stream water is graded as follows: Extremely poor (IWQ ≤ 19), poor (19 < IWQ ≤ 36), regular (36 < IWQ ≤ 51), good (51 < IWQ ≤ 79), excellent (79 < IWQ ≤ 100). It is worth to note that the IWQ index is rather sensitive to small changes in the bearing parameters, given the multiplicative formulation of Equation (1). As corollary of this conception, a good water quality (IWQ > 51) requires that all q values are high, while an excellent quality (IWQ > 79) implies that all q scores are very high.

#### *2.4. Multiple Linear Regression with Interactions*

A typical multiple linear regression model involving a dependent variable (*Y*) and two independent variables (*X*<sup>1</sup> and *X*2) is written as:

$$Y = \beta\_1 X\_1 + \beta\_2 X\_2 \tag{2}$$

where β<sup>1</sup> and β<sup>2</sup> are regression coefficients representing the main effects of *X*<sup>1</sup> and *X*<sup>2</sup> on the predicted values of *Y*. Sometimes, besides the main effects, the dependence of *Y* on *X*<sup>1</sup> and *X*<sup>2</sup> is further constrained by interaction effects. An interaction exists when the effect of one independent variable changes, depending on the value of the other independent variable. In those cases, Equation (2) is recast as [46]:

$$Y = \beta\_1 X\_1 + \beta\_2 X\_2 + \beta\_3 X\_1 X\_2 \tag{3}$$

where β<sup>3</sup> is the joint effect of *X*<sup>1</sup> and *X*<sup>2</sup> on *Y* and the product *X*1*X*<sup>2</sup> is the interaction between *X*<sup>1</sup> and *X*<sup>2</sup> producing that effect. This product is also called a two-way interaction term, because it is the interaction between two independent variables.

The presence of interactions in multiple regression can be identified through statistical tests, namely through assessing the statistical significance of the interaction term, and comparing the coefficient of determination with and without the interaction term. If the interaction term is statistically significant, the interaction term is probably important. And if the coefficient of determination is also much bigger with the interaction term, it is definitely important. If neither of these outcomes is observed, the interaction term can be removed from the regression equation. As alerted in various forums (e.g., https://statisticsbyjim.com/regression/interaction-effects/), when interaction effects are present, it means that interpretation of main effects without considering joint effects may be incomplete or misleading.

In the present study, the independent variables used to model water quality with multiple regression were the riparian buffer strip width and the ratio between forest and agriculture, represented by the variables *L* and *F*/*A* in Table 1. In the first run, the dependent variable was the water quality index represented as IWQ in Table 3. In a second run, the regression analysis was replicated for water temperature, pH, dissolved oxygen, turbidity and total dissolved solids, which are the formation parameters of IWQ, to evaluate their specific roles in the studied area. The dataset comprehended the values of all these variables in the eight catchments. These data were processed for estimation of main effects and joint effects in the *STATISTICA* computer program (http://www.statsoft.com).

#### *2.5. Thematic Maps*

The thematic maps (e.g., Figures 1–3) were prepared in ArcMap software of ESRI [47], a common tool in spatial analysis of hydrologic and environmental data widely used in many studies [42,44,48–66]. The base information was compiled from various spatial databases, namely the maps published by the Brazilian Institute for Geography and Statistics (https://ww2.ibge.gov.br) on the 1: 100,000 scale, and the digital terrain model obtained from the ASTER GDEN V2 satellite image with a spatial resolution of 30 m.

#### **3. Results**

The scatter plots representing the IWQ index as a function of variables L (buffer strip width) and F/A (ratio forest over agriculture), as well as of interaction term L× F/A, are displayed in Figure 6a, 6b and 6c, respectively. The corresponding coefficients of determination are *R*<sup>2</sup> = 0.61, *R*<sup>2</sup> = 0.72 and *R*<sup>2</sup> = 0.93, meaning that the variance explained by the models raises from the main effects (L, F/A) to the interaction effect (L × F/A). The scatter plot of Figure 6a may be influenced by the small number of buffer strip widths (just three different values). This may limit the use of buffer strip width as a continuous variable in a regression model. The results of multiple regression support the conclusions taken from observation of Figure 6. In Table 4, it is evidenced that all multiple regression terms are significant at *p* <sup>≤</sup> 0.05 and that the adjusted coefficient of determination is slightly higher (*R*<sup>2</sup> = 0.99) when the interaction term is incorporated in the model, relative to the no interaction case (*R*<sup>2</sup> = 0.97).

**Figure 6.** IWQ scatter plots. (**a**) This variable is projected as a function of buffer strip width (L), (**b**,**c**) as a function of ratio forest over agriculture (F/A) and interaction term L × F/A, respectively. The data used to draw the plots are depicted in Table 1 (L, F/A and L × F/A) and 3 (IWQ).


**Table 4.** Results of multiple regression model. In this run, the IWQ was used as a dependent variable.

In keeping with the results of multiple regression (Table 4), the relationship between water quality (IWQ), buffer strip width (L) and landscape composition (F/A) can be expressed by the following Equations:

$$INQ = 26.14 + 0.08 \times L + 5.64 \times \frac{F}{A} \,\text{'}\tag{4a}$$

$$IWQ = 23.80 + 0.16 \times L + 10.31 \times \frac{F}{A} - 0.15 \times L \times \frac{F}{A'} \tag{4b}$$

where Equation (4a) represents the relationship without considering the interaction between *L* and *F*/*A* and Equation (4b) the case where this interaction is accounted for. The graphical representation of Equation (4a,b) are illustrated in Figure 7a,b. Figure 7a portrays a couple of parallel lines describing the evolution of *IWQ* as a function of *L*, for the extreme values of *F*/*A* (0 and 1). The lines are parallel because the model predicts no interaction [67]. The buffer strip widths required to ensure a regular water quality in the studied sub-basins (36 < *IWQ* ≤ 51) range from *L* = 130 m to *L* = 310 m, when *F*/*A* = 0, and from L = 60 m to *L* = 250 m when *F*/*A* = 1. These results change considerably when the interaction term is incorporated in the regression analysis, as demonstrated in Figure 7b. Now, the non-parallel lines indicate much smaller *L* ranges to attain the regular water quality, namely *L* = 75–175 m for *F*/*A* = 0, and *L* = 45–155 m for *F*/*A* = 1. The consequences for planning of adopting one or the other model are evident, either for the planning of landscape composition or buffer strip widths. The interaction model may be more realistic because of its larger R2 and interaction term significance.

**Figure 7.** Interaction multiple regression plots: (**a**) No interaction model; (**b**) interaction model.

The results of multiple regression applied to the IWQ parameters are depicted in Table 5. Only the regressions with interaction term were considered in this second run. The results suggest a dominance of turbidity in the control of IWQ in the studied sub-basins. For this parameter, the coefficient of determination (R2 = 0.6) is satisfactory, and all regression coefficients are significant at *p*-level <sup>≤</sup> 0.05. The results for total dissolved solids are characterized by a moderate R<sup>2</sup> = 0.5, but the regression coefficients are not significant. The results for the other parameters are characterized by a low R2 = 0.1 and non-significant regression coefficients.


**Table 5.** Results of multiple regression model. In this run, the formation parameters of IWQ were used as dependent variables, namely water temperature, pH, turbidity, dissolved oxygen and total dissolved solids. The significant values (*p*-level ≤ 0.05) are represented in boldface.

#### **4. Discussion**

The results of multiple regression, with and without interaction effects (Equation (4a,b)), indicate the improvement of water quality as a function of increasing forest to agriculture ratios (landscape composition) and buffer strip widths. The non-scaled regression coefficients (β; Table 4) point to a 45% contribution of L and 55% of F/A to IWQ values in the studied catchments, when the non-interaction model is used, and an equal contribution (50%) from both variables when the interaction model is adopted. These results are not surprising because many studies so far reported the benefits of forest cover and buffer strip widths to stream water quality, as quoted in the Introduction section [8–12].

The striking result refers to the interaction between F/A and L (Equation (4b)), because it describes a substantial reduction of L values required to sustain a regular water quality (36 < IWQ ≤ 51) in the streams, relative to the non-interaction model. The reductions related to the interaction reach 44% when the lines representing F/A = 0 are compared for L values across Figure 7a,b, and 38% in the case of F/A = 1. This is remarkable and requires a thorough interpretation. It seems that the positive effects

on water quality resulting from the independent actions of F/A and L are amplified by a combined action (F/A × L) that ought to identify and justify.

The runoff that reaches a riparian buffer is mostly generated upstream within the catchment. It is, therefore, acceptable to link the combined action to hydrologic processes taken place away from the buffer strips, and hence, related to F/A. These processes must, however, describe a hydrologic connection between forested areas (high F/A) and buffer strips, because a combined action inherently assumes an interplay between the two land cover parameters. A potential strong candidate is the infiltration capacity of soils, and, more importantly, its spatial distribution [68]. This is a crucial control on hydrological processes at the interface between the ground surface and soil, including the distribution of flowing water by overland flow and shallow underground flow. The continuity of overland flow or alternation with subsurface flow depends on the intensity of rainfall and the locations of relatively high infiltration patches isolated on hillslopes [69–71]. In this context, it is expected that high F/A catchments comprise a larger number of infiltration patches, and that subsurface flow dominates in these catchments. It is also expected that subsurface flow is widespread within the soil layer, a condition that improves the probability of soil water to cross a buffer strip before reaching the stream. The higher the F/A ratios, the larger will be the chance of soil water to interact with the buffer strip. This scenario could explain a combined action of F/A and L on IWQ. Contrarily to high F/A catchments, low F/A catchments (e.g., dominated by agriculture) are prone to the generation of an overland flow network, because the absence of permanent vegetation substantially reduces the number of infiltration patches. The overland flow network channelizes runoff and conveys the surface water into specific confluence points within the stream, reducing or even hampering an interaction with a buffer strip [72]. Taken altogether, the F/A and L parameters represent the buffering capacity of vegetated areas distributed away and along the stream banks, respectively, while the F/A × L parameter represents the additional capacity promoted by a hydraulic connection between the two areas.

Besides the infiltration capacity issue, it is worth to explore the potential influence of rainfall intensity in the interaction effect, considered relevant in areas where annual climate fluctuations are stronger, including the tropical regions [73,74]. Turbidity is extremely sensitive to rainfall intensity and played a dominant role in the regression analyses of individual parameters (Table 5). The turbidity records of all studied catchments (see Supplementary Materials) are represented graphically in Figure 8, along with the corresponding daily precipitation record. It is evident the response of turbidity to larger values of daily precipitation, namely in the Borá 1 (January 2016), Lageado (April 2016) and Alegria (January 2017) catchments. A close inspection of Table 1 reveals that these catchments are among those with a lower L × F/A value. There are, however, striking exceptions: The Uberaba River catchment has the lowest L × F/A value (1.3 m), but barely responds to precipitation events; the Borá 2 catchment turbidity peak in September 2016 is not justified by a precipitation event. It is worth recalling, however, that the control of turbidity in streams is not only determined by the studied parameters, namely F/A, L, infiltration capacity and flow convergence related to L × F/A, and rainfall intensity. Topography is also a key parameter [75], which was not addressed in this study because the focus here was put on the interaction between landscape composition and buffer strip widths. We believe topography would help to explain the Uberaba River exception. This is a headwater catchment located in a relatively flat area (compare Figure 4 with Figure 1). These areas are well acquainted with retain sediments because they generate lower overland flow velocity while maximizing infiltration and particle deposition, in opposition to steep slope areas [76].

The regression results based on individual parameters (Table 5) exposed a significant relationship between catchment variables (landscape composition, buffer strip width), including their interactions, and water turbidity, but did not reveal identical influences of those variables or their interactions on other parameters measured in water. It should be remembered, however, that water quality parameters may respond differently to catchment variables depending on the spatial scale or antecedent rainfall conditions, as noted in Uriarte [31]. We, therefore, clarify that our results are valid at the studied

spatial scales (Figure 4) and antecedent rainfall conditions (Table 2). Transposition to other settings needs verification.

**Figure 8.** Values of daily precipitation and stream water turbidity observed in the studied catchments during the sampling campaigns. The discrete values are provided as Supplementary Material.

The recognition of an interaction effect between landscape composition and buffer strip widths capable of amplifying stream water quality improvements is certainly beneficial for water resources planning and management. The dominant role of turbidity in the regression analyses of individual terms suggests that water quality deterioration in the studied sub-basins is mostly related to soil erosion and sediment transport rather than with leaching of dissolved fertilizers from the arable land into the stream. Therefore, conservation measures related to the control of soil erosion should be prioritized in this protected area. Eventually, the most cost-effective measure is the reinforcement of manuring and composting of soil to raise the levels of organic matter [77] and produce stable aggregates that are resistant to detachment by rain drop action. The level plantation is another measure of soil erosion control, frequently used to reduce laminar erosion and increase soil water uptake. The level plantation can be successfully coupled with strip cropping [78] that involves the alternation of forages with strips of row crops, such as sugar cane. The control of soil erosion also comprises implementation of techniques that prevent soil compaction, such as no-tillage treatments [79], maintenance of crop residue to keep organic matter and nutrients in the topsoil, or more costly measures, such as terraces in level or gradient, since they reduce the ramp length reducing the surface drag of particles and nutrients [80]. To become effective, implementation of conservation measures should be monitored within spatial decision support systems focused on water resources planning and management [81], and integrate public policies and environmental management plans.

The regression results showed that, even considering the interaction effect regular water quality in the studied catchments is only attained when buffer strip widths are within the ranges 75 ≤ L ≤ 175 m (F/A = 0) or 45 ≤ L ≤ 155 m (F/A = 1) (Figure 7). These values are much larger than legal values in force in Brazil. So far, the Brazilian law has defined the buffer width limits based on two scenarios: The Old Forest Code (Federal Law No. 4771/1965) and the New Forest Code (Federal Law No. 12651/2012). In the first case, for watercourses up to 10 m wide the permanent preservation areas needed to extend at least 30 m upwards from the stream margin considering the widest seasonal riverbank. In the second case, there are two rules: The transition rule takes into account the size of land property calculated as fiscal modules and creates a distance from the stream margin that goes from a minimum of 5 m to a maximum of 20 m, considering the regular river bank; the permanent rule defines 30 m as unique distance. This study reinforces the suggestion of Valera et al. [15], who alerted that a 30 m buffer strip width, as proposed in the New Forest Code, is barely capable of protecting water quality in the EPA-MURB. The discussion on buffer strips, their geometry and composition, optimal

widths, cost-benefit analysis for implementation [82–85], among other topics, is still a matter of debate. The discussion on interaction effects is expected to become another topic on this so challenging analysis.

#### **5. Conclusions**

The results of a multiple regression model involving an interaction term revealed the combined positive influence of landscape composition and buffer strip widths (L) on stream water quality, in eight catchments of Uberaba River Basin Environmental Protection Area (Minas Gerais State, Brazil). The landscape composition was characterized by the forest to agriculture ratio (F/A), and high F/A catchments were viewed as basins where forested areas located away from the stream are in hydraulic connection with riparian vegetation distributed along the stream banks. This hydraulic connection presumably amplifies the buffering capacity of forested areas and riparian buffers acting independently. To our best knowledge, this is a new finding in the study of buffer strips and their relationship with stream water quality. In practice, the combined effect reduced the width of buffer strips required to keep water quality at a fair level. Without the interaction, the calculated L range was 60–310 m. It decreased to 45–175 m when the interaction term was accounted for in the regression model. The reduction was expressive, but not enough to drop below the maximum legal value imposed by the Brazilian Forest Code (30 m). Eventually, this legal framework could be adapted to our scientific findings. The problem of setting thresholds to buffer strip widths is not exclusive from Brazil, and therefore, our results could serve as alert to other national realities.

#### **Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/9/1757/s1.

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

**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); Agência do Ministério da Ciência, Tecnologia, Inovações e Comunicações (MCTIC); and Land Use Policy Brazilian Group (POLUS). The author affiliated to IFTM Renato Farias do Valle Júnior wishes to acknowledge the funding through the CNPq research scholarship Proc. 307921/2018-2. For the author integrated into the CITAB research center, it was further financed by the FEDER/COMPETE/POCI—Operational Competitiveness and Internationalization Program, under Project POCI-01-0145-FEDER-006958, and by National Funds of FCT—Portuguese Foundation for Science and Technology, under the project UID/AGR/04033/2019. For the author integrated into the CQVR, the research was additionally supported by National Funds of FCT–Portuguese Foundation for Science and Technology, under the project UID/QUI/00616/2019.

**Acknowledgments:** Mauro Ferreira Machado is acknowledged for fruitful discussions and sharing of information.

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

#### **References**


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

## *Article* **A Methodology for the Fast Comparison of Streamwater Diurnal Cycles at Two Monitoring Points**

#### **Andrei-Emil Briciu 1,\*, Adrian Graur 2, Dinu Iulian Oprea <sup>1</sup> and Constantin Filote <sup>2</sup>**


Received: 23 October 2019; Accepted: 27 November 2019; Published: 29 November 2019

**Abstract:** There are numerous streamwater parameters that exhibit a diurnal cycle. However, the shape of this cycle has a huge variation from one parameter to another and from one monitoring point to another on the same river. Important variations also occur at the same point during some events, such as high waters. Water level, specific conductivity, dissolved oxygen, oxidation reduction potential, and pH of the Suceava River were monitored for 365 days (2018–2019, hourly sampling frequency) in order to assess the upstream-downstream changes in the diurnal cycle of these parameters, some of these changes being caused by the impact of Suceava city, which is located between the selected monitoring points. The multiresolution analysis of the maximal overlap discrete wavelet transform and the wavelet coherence analysis were combined in a flexible methodology that helped in comparing the upstream and downstream shapes of the diurnal cycle. The methodology allowed for a fast comparison of diurnal profiles during periods of high waters or baseflow. Notable changes were observed in the moments of diurnal maxima and minima.

**Keywords:** script files; mean diurnal profile; wavelet coherence

#### **1. Introduction**

The wavelet analysis of streamwater parameters has become more and more popular in the last two decades due to the advantages of this method for the study of non-linear processes [1,2]. Wavelet analysis techniques were applied mostly on stage or discharge time series due to the wide availability of this type of data [3]. It is in the last decade when water quality parameters were involved intensively in wavelet analyses [4,5]. Some physical and chemical properties of streamwaters are proper indicators of water quality and can be used to trace the environmental impact of man. The specific conductivity, dissolved oxygen, oxidation reduction potential (ORP), and pH of a river water can be used to indicate the impact of urban areas on the environment and are sometimes included in wavelet analyses [5–7]. Cities modify the properties of natural waters through various active or passive processes, such as the discharge of stormwater runoff or the creation of an urban heat island, which has multiple effects on urban waters [8]. Urban wastewater treatment plants can alter the diurnal profile of a streamwater chemistry parameter [9,10].

The diurnal cycle in streamwaters is caused by the Earth's rotation around its axis, which leads to the day/night cycle. This cycle modifies evaporation and evapotranspiration in catchments [11,12]. These changes can be measured as variations in the temporal evolution of numerous streamwater parameters, which often have an interdependent behavior with the diurnal oscillation [13]. Natural or anthropogenic events (such as rainfall or pollution events) add transient fluctuations in the diurnal cycle. Finding a relevant shape of the diurnal cycle is needed in order to distinguish the periodic and aperiodic modifications in case study time series. A high resolution diurnal profile is obtained from

high frequency measurements. This is why, up-to-date, such profiles are missing in Romania for most water quality parameters. Existing studies on diurnal streamwater profiles have focused on water level or water temperature [10,14].

Streamwater monitoring in Romania is done mainly by the Romanian Waters National Administration, whose data is sometimes used in studies about water chemistry of the Suceava River [15], but this data lacks high frequency sampling and the needed spatial density along a river. The environmental impact of some water contaminants is often assessed for various areas in Romania [16–20] and cities have a traceable impact on streamwater quality [21].

A monitoring system that measures Suceava River water quality upstream and downstream of Suceava city was implemented in 2018 and this study aimed to provide a methodology for a fast analysis of the diurnal cycles of various streamwater parameters at two monitoring points. The wavelet analysis was used to reveal the spatial and temporal variations in the streamwater diurnal cycles. The speedier analysis is needed for the growing sizes of databases. It is often useful to apply an analysis method that acts as an easily adjustable preview of data in order to identify interesting phenomena for further analyses. Also, to our knowledge, this was the first study that computes the average diurnal profiles of specific conductivity, dissolved oxygen, ORP, and pH for the Suceava River.

The methodology proposed in this study was applied to streamwater data from monitoring points located above and below a city in order to discover some details of how the diurnal profile of a river is affected by a city. In this paper, we indicate some links between the observed variations in diurnal profiles and urban elements, such as the urban wastewaters or the urban heat island, that generate the environmental impact of a city.

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

#### *2.1. Study Area*

The monitoring system. which provided data in this study, was implemented by the University of Suceava and consisted of two monitoring points where water properties are measured every hour (the monitoring points are 11.6 km in a straight line away from each other, with Suceava city in the middle of this distance). Suceava city has a total number of inhabitants fluctuating around 100,000 people.

The study area was located in Suceava Plateau, part of the Moldavian Plateau. Climate is temperate continental with warm and wet summers (the average annual air temperature is 7.86 ◦C and the annual sum of precipitation is 578 mm [10]). The streambed elevation between the selected points ranged approximately from 280 to 260 m above sea level (a.s.l.) and has a sinuous path, which helps in mixing streamwater so that it records identical values at any point of a transversal section. The upstream monitoring point was at Mihoveni Dam (47.681◦ N, 26.2◦ E). This is a quasi-inactive dam aimed at generating run-of-the-river hydroelectricity; it operates only during some high water events for regulating the water level. The downstream monitoring point was placed at Tis, ăut,i (47.618◦ N, 26.323◦ E), downstream of an urban wastewater treatment plant and some floodplain landfills (Figure 1). The Suceava River has an average flow rate of 16.87 m3/s inside Suceava City [10]. The Suceava River has two periods of high waters. The first one occurs during the middle of springtime (April, monthly average flow rate of 29.7 m3/s) due to snowmelt in the catchment, especially in the mountain area, while the second one happens at the beginning of summer (June, 29.9 m3/s) as result of heavy rainfalls [10]. During high water events, the Suceava River frequently exceeds 50 m3/s. The lowest monthly average flow rate is recorded during the winter (January, 6.2 m3/s) because of the negative air temperatures and little precipitation [10]. The urban tributaries of the Suceava River in the study area had a total discharge of approximately 0.5 m3/s. The water chemistry of rivers in north-eastern Romania, which includes our study area, was briefly discussed in some studies [22,23], some of them highlighting the intense self-purification processes of Suceava River water inside and downstream of the city [24]. However, these studies were not based on long time series with high frequency of measurements.

**Figure 1.** Map of the study area, including the location of the monitoring points and the urban wastewater treatment plant (WTP).

#### *2.2. Instruments*

Two AquaTROLL 500 multiparameter probes were used to monitor physical and chemical streamwater properties (paired with two Tube 300R for telemetry). The instruments were equipped with sensors that measured: pressure/level (accuracy: ±0.1% full scale (9 m), resolution: 0.01% full scale), electrical conductivity (automatically converted to specific conductivity (SC); accuracy: ±0.5% of reading +1 μS/cm, resolution: 0.1 μS/cm), dissolved oxygen (DO; accuracy: ±0.1 mg/L, resolution: 0.01 mg/L), oxidation reduction potential (ORP; accuracy: ±5 mV, resolution: 0.1 mV), and pH (accuracy: ± 0.1 pH unit, resolution: 0.01 pH).

#### *2.3. Data*

Data analyzed in this study were recorded from 10 October 2018 until 9 October 2019 (with one exception: DO data is missing from November 7 to December 5, 2018). The local standard time (EET) was used for all time series. Corrections were applied to our data as follows: compensations were applied after periodic instrument calibrations, missing singular values were obtained by using linear interpolation, outliers were removed through replacement with values from data smoothed with a moving average filter, and some missing consecutive values were obtained by using regression equations. At the downstream sampling point, the measurements were sometimes affected by residuals from wastewaters and the time series were verified by using independent measurements carried out using a Hach HQ40d portable multiparameter instrument. Data of all parameters can be viewed and downloaded at the project website, http://water.usv.ro/data.php (where a map of the study area can also be inspected).

The yearly average of the specific conductivity was higher downstream (at Tis, ăut,i; 549 μS/cm) than upstream (at Mihoveni; 483.1 μS/cm) and the standard deviation values had the same relationship (97.9 μS/cm versus 72.4 μS/cm). Urban wastewaters were the main cause of the differences; treated and untreated waters are discharged into the Suceava River through the wastewater treatment plant effluent and Cetăt,ii Creek [10,23,25]. The difference is much larger after snowfalls, when de-icing measures are applied and/or when high air temperatures lead to important snowmelt and runoff from

roads and roofs (see, for example, peaks in SC that were only recorded in the second half of November 2018 downstream after large snowfalls occurred) (Figures 2 and 3). Dissolved oxygen had mean values of 10.6 and 8.82 mg/L (upstream and downstream, with corresponding standard deviations of 2.06 and 2.69 mg/L, respectively). The lower average value downstream of the urban areas was an effect of a warmer and more polluted streamwater. The urban heat island was observed in the study area in a previous study [10].

ORP mean values upstream and downstream were 412.16 and 338.01 mV and standard deviation at these points was 32.11 and 52.97 mV, respectively. The values of pH were 8.45 (upstream) and 8.22 (downstream) and the difference in standard deviations was the greatest of all parameters: 0.09 (upstream) and 0.31 (downstream).

**Figure 2.** Variations in the studied parameters of the Suceava River at the upstream monitoring point (Mihoveni) from 10 October 2018 to 9 October 2019: (**a**) specific conductivity (SC), (**b**) dissolved oxygen (DO), (**c**) oxidation reduction potential (ORP), and (**d**) pH.

**Figure 3.** Variations in the studied parameters of the Suceava River at the downstream monitoring point (Mihoveni) from 10 October 2018 to 9 October 2019: (**a**) specific conductivity (SC), (**b**) dissolved oxygen (DO), (**c**) oxidation reduction potential (ORP), and (**d**) pH.

The average diurnal profiles of SC, DO, ORP, and pH indicated differences in the shapes and positions of the diurnal maxima and minima not only between parameters, but also between the values of the same parameter at the two monitoring points (Figure 4). SC upstream had a diurnal maximum during midday (11:00–12:00), while downstream, the maximum was recorded at 03:00. The maximum values of DO recorded were in the late afternoon and early evening at the selected monitoring points; this difference was not caused by the streamwater temperature, which had similar moments of maxima in the study area [10]. Rather, this similarity happens instead of an inverse relationship that should occur theoretically. The average hourly values of downstream ORP created a prolonged interval with maximum values in the first third of the day. The diurnal minimum of upstream pH occurred at 13:00 and might be partly explained by the theoretical inverse relationship between pH and temperature, but the downstream pH maxima was recorded at 15:00.

The average diurnal profile of streamwater level at both monitoring points had irregular fluctuations that do not satisfy the definition of a diurnal cycle. This was due mainly to rainfall events, which have strong responses in the change of the river water level. If the analysis window was restricted only to periods of baseflow, diurnal cycles could be observed for this parameter too (minima in the midday, maxima during late evening). The mean relative streamwater level was

0.786 m upstream of Suceava city and 0.895 m downstream. The standard deviation indicated a higher difference between the two points: 0.121 m upstream and 0.305 m downstream. The increased variability of the water level downstream of Suceava city is to be attributed mainly to water runoff from impervious urban areas because of tributaries that flow into the Suceava River between the two sampling points. These tributaries have small discharges at baseflow and receive some of the pluvial drainage from the metropolitan area (especially S, cheia, Cetăt,ii, Dragomirna, and Podu Vatafului creeks [9]).

**Figure 4.** Average diurnal profiles of the selected parameters at the upstream (1) and downstream (2) monitoring points (horizontal axes in hours): (**a**) specific conductivity, (**b**) dissolved oxygen, (**c**) oxidation reduction potential, and (**d**) pH (data from 10 October 2018 until 9 October 2019; with the exception of DO, which was calculated from 6 December 2018 until 9 October 2019 due to missing data).

#### *2.4. Analysis Methods*

Variations in flow rates/water level can cause important changes in the diurnal profile of other streamwater parameters. Increases in water level have a clear cause, rainfall, and the slow decline in water level is to be naturally expected until the next rainfall. For this reason, we chose the water level parameter as the control factor when analyzing SC, DO, ORP, and pH. More precisely, we analyzed the water level time series in order to choose the position and length of the analysis windows used for comparing the diurnal profiles of the other parameters. Because the increase in water level is the active change, we decided to split the water level time series into subsets whose ends were given by the minimum values between high waters (Figure 5). A time interval between such minima would be a good case study to observe changes in diurnal cycles induced by peaks in flow rates.

**Figure 5.** Suceava River water level at the (**a**) upstream and (**b**) downstream monitoring points plotted by using (1) raw data and (2) simplified data computed with the multiresolution analysis. Vertical red lines represent the limits of time intervals between relevant minima (10 October 2018 –9 October 2019).

Data analysis consisted of four major steps: (1) applying the multiresolution and wavelet analysis on raw level data in order to obtain time intervals for case studies; (2) smoothing the time series of water level, SC, DO, ORP, and pH, obtaining the evolution of diurnal cycles with the super-daily trends removed and normalizing the remnant time series for the upstream-downstream comparison; (3) computing the average diurnal profiles for case study time intervals from the normalized time series; and (4) comparing the case study time intervals from the normalized time series through the wavelet coherence analysis.

Prerequisites for the proposed analysis methodology were: MATLAB software (for all steps), the three script files found in the Supplementary Material of this paper (for the first three major steps), the cross wavelet and wavelet coherence toolbox of Aslak Grinsted [26] (https://www.mathworks. com/matlabcentral/fileexchange/47985-cross-wavelet-and-wavelet-coherence), data loaded/created in MATLAB workspace consisting of two variables named "upstream" and "downstream", which were 1D column vectors representing one parameter per monitoring point and having a length equal to a multiple of 24 (because data in this study was sampled 24 times per day).

The first step had four stages:

• Applying the maximal overlap discrete wavelet transform to the water level data from both monitoring points. This was acquired with the function *modwt*, which has the following syntax:

$$y = 
mod wrt(x, 
wname, n)\_r$$

where *x* is the real-valued signal; *wname* is the name of an orthogonal wavelet, which, in our case, was '*haar*'; and *n* represents the desired number of levels of detail coefficients—7 in this case;

• Executing the multiresolution analysis based on *modwt*, with the function *modwtmra*, which had a similar syntax:

$$z = 
mod w 
times 
mu(y, 
w 
mu 
me)\_r$$

where *wname* should be the same as that used with *modwt*;


The wavelet decomposition and the maximal overlap discrete wavelet transform were applied previously in a few studies concerning some streamwater parameters, such as water level, DO, and pH [3,5,27,28], which also offer the mathematical description of the mentioned wavelet analyses. Haar wavelet was successfully tested for water parameters [29].

The second step had 3 stages:

• Smoothing the raw time series by using the smooth function, which had the following syntax:

#### *yy* = *smooth(x,span,method)*,

where the *method* used here was the moving average (*moving*; lowpass filter with coefficients equal to the reciprocal of the span) and the *span* was equal to a day (24). Note that, due to software limitations, span was reduced by 1 to an even number;


Steps 1 and 2 were computed by executing the *analyze* command line in MATLAB console. This line applies the codes written in the analyze.m script file found in the Supplementary Material. After this command line was applied, two new variables, normalisedU and normalisedD, were created in the workspace and represented the normalized data from upstream and downstream, respectively, obtained at step 2, stage 3 (some other intermediary variables were also created in the workspace, but were handled automatically by codes). The *analyze* command line also generated a graphical output, similar to Figure 5, and was applied to water level data. In order to obtain only normalized data (as in step 2), for the other streamwater parameters, the command line *transform* had to be executed, which used the transform.m script file (in the Supplementary Material) and generated the variables PnormalisedU and PnormalisedD.

All of the new time series from (P)normalisedU and (P)normalisedD were involved in step 3, because they served as the source of data that was cut for case studies by using the time intervals established at the end of step 1. Data selected for case studies had to be included in two new variables, named CSnormalisedU and CSnormalisedD.

Step 3 had only one stage, when the average diurnal profile was computed, and some graphical outputs were generated, as in Figure 6 (using the newly generated variables CSprofileU and CSprofileD). This step was executed with the command line *compare* (with codes written in the compare.m script file, found in the Supplementary Material).

Step 4 had one stage and used the *wtc* function which performs the wavelet coherence analysis between two variables and has the syntax:

#### *wtc(a,b),*

where *a* and *b* are values from the upstream and downstream parameters, respectively (e.g., normalisedU, PnormalisedU, or CSnormalisedU). The analysis produced scalograms that indicated how strong two

signals co-vary by using a power spectrum and phase arrows. Also, periodicities in signals were tested against the AR1 red noise through a Monte Carlo test and the relevant results were marked with a black line. The wavelet coherence used the Morlet mother wavelet because it is best suited for comparing time series with non-linear processes [26,30]. Wavelet coherence mechanics and results were described by previous studies that applied this type of analysis to hydrological time series [26,31]. The script files used in this methodology are included in the Supplementary Material and described in an accessible manner between the lines of codes contained in the files.

**Figure 6.** Suceava River water level at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (17 June–9 July 2019).

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

The analysis of the major convolutions/waves of water level during a year led us to select two subsets of time series that should serve as relevant case studies of changes in diurnal profiles depending on water level variability. Case study 1 comprised 23 days during a period of intense rainfalls and high waters (17 June–9 July 2019; Figures 6–10), while case study 2 comprised 10 days during the autumn baseflow (20–29 October 2018; Figures 11–15). The length of the subsets with common major evolutions could be increased or decreased depending on the number of levels of detail coefficients of the maximal overlap discrete wavelet transform, which in our case was set to 7.

Depending on the characteristics of water flow or the need for longer/shorter subsets, the number of levels of detail coefficients could be edited in the anlyze.m script file. The smaller the number of levels, the fewer details that were removed from the original time series through decomposition. The subsets were very numerous and small; more levels removed more details and generated and abstracted the approximation of the original data, where only a few subsets could be assessed. Therefore, the choice of decomposition level depends on the user/author [32], as it is the choice of applying the multiresolution analysis only on the water level or on other parameter/all parameters. Available orthogonal wavelets in MATLAB that can replace the Haar wavelet (*'haar'*) option are: Daubechies (*'dbN'*), Symlets ('*sym4*' or *'symN'*), Coiflets (*'coifN'*), and Fejér–Korovkin (*'fkN'*).

**Figure 7.** Suceava River SC at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (17 June 17–9 July 2019).

**Figure 8.** Suceava River DO at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (17 June 17–9 July 2019).

Distinct average diurnal profiles can be observed for any water parameter in case study 1, except for SC (Figure 7), which was very sensitive to water level oscillations during high waters. Such time intervals explain the complex shape of the average diurnal profile of SC calculated for a year (Figure 4a1). Another parameter that had the diurnal cycle easily altered by important changes in water level was pH (Figure 10). Changes in the shapes and positions of the diurnal cycles occurred in the average diurnal profile of DO (Figure 8); at both monitoring points, diurnal maxima moved towards midday or late morning (Figure 4b).

**Figure 9.** Suceava River ORP at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (17 June 17–9 July 2019).

**Figure 10.** Suceava River pH at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (17 June 17–9 July 2019).

The script files analyze.m and transform.m allowed for editing the smoothing and detrending methods. The smoothing can be done with a filter other than the default moving average. Available filters were: '*lowess*', a local regression with weighted linear least squares and a 1st degree polynomial model; '*loess*', the same as the previous, but with a 2nd degree polynomial model; '*rlowess*' or '*rloess*', robust versions of '*lowess*' or '*loess*' with lower weight for outliers in the regression; '*sgolay*', a Savitzky–Golay filter. Also, smoothing could be set to be more or less aggressive with the increase or decrease in the span number.

The detrending technique may be different than the proposed difference between the raw and smoothed date. A method with a similar ability to remove general and seasonal trends might show the difference between adjacent elements, if the following syntax is used: *y*=*di*ff*(x)*.

We inserted a few codes between the lines that execute smoothing and detrending. These are aimed at helping those that intend to obtain detrended values that can no longer be normalized. Instead, a detrended time series with absolute values, similar to those of the raw time series (but lower), will be obtained. It is necessary that this option is applied only to data with inter-diurnal variations significantly higher than the annual oscillation and with weak variability in the diurnal cycle (dedicated variables created in the workspace were differenceX and subtrendX, where X is U (upstream data) or D (downstream data)).

The time interval of case study 2 belonged to the autumn baseflow, where the main aperiodic changes were caused by a few rainfalls. For the first time, SC average diurnal profiles at both monitoring points were similar (Figure 12). This led to the conclusion that, during periods of baseflow, which are longer than the periods affected by intense rainfalls, this similarity is frequent and only these average diurnal profiles might be considered relevant for describing various persistent riverine processes.

Another interesting observation was linked to the diurnal behavior of DO at the downstream monitoring point (Figure 13b). The diurnal cycle had the same relative position of minima as observed in the annual diurnal profile (Figure 4b2), but recorded maxima in two secondary peaks superposed on the diurnal oscillation. Various processes may co-generate this double-peaked shape, including the urban heat island (because the water is measured at the downstream point after passing through the middle of Suceava city). A change also occurred in the hourly position of the maximum value of the DO average diurnal profile: it occurred early in the morning, instead of late evening, as in the annual profile.

**Figure 11.** Suceava River water level at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (20–29 October 2018).

The diurnal oscillations of ORP and pH were strongly affected by rainfalls during baseflow at the upstream monitoring point. This was caused by the weak intensity of the diurnal cycle, which became easily modifiable by quasi-random events. In contrast, at Tis, ăut,i, the diurnal cycles of ORP and pH had oscillations stronger than the general trend, which was different at the downstream point than at the upstream point. Given the fact, discussed previously, that the metropolitan area modifies the averages and standard deviations of the measured parameters, we may assume that the stronger diurnal cycle of ORP and pH at the downstream point is regulated by Suceava city. This may also explain why pH values were not lowered by rainfalls, which, by contrast, cause relatively sudden drops at the upstream point. This latter behavior is natural because of the acidic nature of the raindrops, which dissolve gases, such as CO2, from the atmosphere. The amplitude of the diurnal cycle of pH had amplitudes much

greater than 0.2 units, which might induce stress to aquatic fauna. The high variability of pH combined with the low DO during some periods of the year could have caused fish mortality, as reported in June 2018 near the downstream monitoring point (at Lisaura). During the period investigated by us in 2018 and 2019, the minimum DO value was 5.2 mg/L upstream of Suceava city and 2.1 mg/L (critical value) downstream of the city. Pollutants from the wastewater treatment plant of Suceava city contributed to the observed values of the monitored parameters, especially when the plant has difficulties in treating the wastewaters according to specific regulations [10]. The changes in the metropolitan area toward higher imperviousness [10] and the changes in climate toward an increased average air temperature in the study area [33] will generate higher pH oscillations and lower DO concentrations downstream of Suceava city and more fish mortality events are to be expected in the near future.

**Figure 12.** Suceava River SC at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (20–29 October 2018)..

**Figure 13.** Suceava River DO at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (20–29 October 2018).

**Figure 14.** Suceava River ORP at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (20–29 October 2018).

**Figure 15.** Suceava River pH at the (**a**) upstream and (**b**) downstream monitoring points plotted as (1) raw data, (2) detrended and normalized data, and (3) the average diurnal profile of the normalized data (20–29 October 2018).

The wavelet coherence analysis allowed us to find time intervals when the diurnal cycle was statistically significant (0.95) at both monitoring points at the same time and to observe the temporal changes of the phasing of this cycle (Figures 16–18; arrows pointing right indicate in-phase evolution and anti-phase when pointing left). At a first glance, we observed the high extension of the area with high power (red) on SC scalograms, very similar to that of the water level scalograms. Except for the diurnal band observed for SC in October 2018 (Figure 17a), the high power indicates only very good co-variance over periods greater than 32/64 h; the weak coherence at high frequencies is caused by the disturbing effect of high waters.

**Figure 16.** Scalograms of the wavelet coherence analysis of some water parameters measured upstream and downstream during 17 June–9 July 2019: (**a**) SC, (**b**) DO, (**c**) ORP, and (**d**) pH.

**Figure 17.** Scalograms of the wavelet coherence analysis of some water parameters measured upstream and downstream during 20–29 October 2018: (**a**) SC, (**b**) DO, (**c**) ORP, and (**d**) pH.

**Figure 18.** Scalograms of the wavelet coherence analysis of the upstream and downstream water level during: (**a**) 17 June 17–9 July 2019 and (**b**) 20–29 October 2018.

The only persistent coherence between time series recorded upstream and downstream in both case studies was observed for DO (Figures 16b and 17b). During June–July 2019, the diurnal cycles at both monitoring points were in phase or almost in phase, while during October 2018 they were in obvious anti-phase on both sides of the cone of influence (which separates time intervals affected by the edge effects from those that are not affected). Moments when other parameters had coherent evolution of the diurnal cycles in both monitoring points were at the end of June and beginning of July 2019 for ORP and in October 2018 for pH.

Wavelet coherence analyses are not useful only for studying the relationship between two points, but also for investigating the inter-parametric links at the same monitoring point [7]. Climate indices and water quality parameters can also be used for wavelet coherence analyses [34].

#### **4. Conclusions**

This study is the first one to use high-frequency (hourly) measurements done over a 365-day interval in the analysis of the temporal and spatial evolution of SC, DO, ORP, and pH of the Suceava River. Raw and normalized time series from case studies, together with average diurnal profiles and scalograms of the wavelet coherence analysis, indicate distinct diurnal evolutions of the selected parameters over time, which cannot be deduced from the annual average diurnal profiles. Important differences in the shape and hourly position of the diurnal cycles of the same parameter occurred between the upstream and downstream monitoring points, despite of the short distance that separated them. Important changes in the monitored water parameters are caused by Suceava city, located between the selected monitoring points.

The case studies time intervals were selected objectively by using the multiresolution analysis of the maximal overlap discrete wavelet transform applied to water level data. The removal of detail coefficients from the raw water level data allowed for detecting the major variations that impose changes in other water parameters. These major variations were generated by increases in water levels due to rainfalls. The time series of SC, DO, ORP, and pH were detrended prior to comparing their average diurnal profiles at the upstream and downstream monitoring points for a better understanding of the diurnal cycles.

The pair of monitoring points, located upstream and downstream of Suceava city, allowed for measuring the impact of urban wastewaters on the Suceava River. The specific conductivity was higher downstream of the city (yearly averages: 483.1 μS/cm upstream, 549 μS/cm downstream) because of both treated and untreated waters discharged directly or indirectly into the Suceava River. Dissolved oxygen mean annual values were 10.6 mg/L (upstream) and 8.82 mg/L (downstream). The low downstream values are caused by pollutants in water, which consume oxygen in chemical reactions, and because of the urban heat island, which led to warmer waters, less able to dissolve oxygen from the atmosphere.

Three MATLAB script files are provided for a fast comparison process, especially useful when using big data. Codes within the script files are explained and possible variations of the analyses are proposed for a more flexible approach when different data are analyzed. Options include, but are not limited to, variations of the wavelet type, decomposition level, and detrending methods. The proposed methodology is best suited for comparing results from case studies.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/12/2524/s1. Supplementary Material (a zip archive containing three MATLAB .m script files).

**Author Contributions:** Conceptualization, A.-E.B. and A.G.; Data curation, A.-E.B. and C.F.; Investigation, A.-E.B. and D.I.O.; Methodology, A.-E.B. and C.F.; Resources, A.-E. B., A.G., and D.I.O.; Supervision, A.G.; Visualization, D.I.O. and C.F.; Writing—original draft, A.-E.B.; Writing—review & editing, A.G.

**Funding:** This research was funded by CNCS - UEFISCDI, project number PN-III-P1-1.1PD-2016-2106.

**Acknowledgments:** Wavelet coherence software was provided by A. Grinsted. This study was supported by data obtained within the research project SQRTDA (Streamwater Quality Real-Time Data Analysis). This work was supported by a grant of Ministry of Research and Innovation, CNCS - UEFISCDI, project number PN-III-P1-1.1PD-2016-2106, within PNCDI III.

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

#### **References**


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

#### *Article*

## **An Improved Model for the Evaluation of Groundwater Recharge Based on the Concept of Conservative Use Potential: A Study in the River Pandeiros Watershed, Minas Gerais, Brazil**

**Marcelo Alvares Tenenwurcel 1, Maíse Soares de Moura 1, Adriana Monteiro da Costa 1, Paula Karen Mota 1, João Hebert Moreira Viana 2, Luís Filipe Sanches Fernandes <sup>3</sup> and Fernando António Leal Pacheco 4,\***


Received: 14 February 2020; Accepted: 28 March 2020; Published: 1 April 2020

**Abstract:** Water resources have been increasingly impacted due to the growth of water demand associated with environmental degradation. In this context, the mapping of groundwater recharge potential has become attractive to water managers as it can be used to direct public policies and conserve this natural asset. The present study modifies (improves) a spatially explicit model to determine groundwater recharge potential at the catchment scale, testing it in the Pandeiros River basin located in the state of Minas Gerais, Brazil. The model is generally based on the water balance approach and the input variables were compiled from institutional sources and processed in a Geographic Information System. The novelty brought by the aforementioned modification relates to the coupling of physical variables (conventional way) and land management practices (introduced here) in the estimation of a percolation factor. The role of land management practices for percolation was assessed by the so-called Conservative Use Potential (PUC) method, which classifies the areas of a river basin in terms of their potential for sustainable use. The results were validated by an independent method, namely the recession curve method based on the interpretation of hydrographs. In general, the groundwater recharge potential is favored in flat to gently undulating areas and forested regions, as well as where the landscape is characterized by well-structured soils, good drainage conditions and large hydraulic conductivity. The map of groundwater recharge potential produced in this study can be used by planners and decision makers in the Pandeiros River basin as a tool to achieve sustainable use of groundwater resources and the protection of recharge areas.

**Keywords:** groundwater recharge; water resources management; Conservative Use Potential; river basin; geographic information system; water balance

#### **1. Introduction**

Clean water resources are becoming scarcer due to the increasing demand for its use and to environmental degradation [1]. It has been estimated that 1/3 of all countries will have to adapt

their productive processes by 2025 because water is lacking, and more than 3 billion people might be living in regions under chronic draughts or water stress [2–6]. Owing to its great importance for society, economy and the environment, it is necessary to properly manage water resources. The increasing degradation of surface water resources not only in Brazil but all over the world [7], makes the sustainable management of groundwater become an even more essential activity. In this context, a precise evaluation of groundwater distribution entails the understanding of aquifer refilling processes and the quantification of groundwater inputs. This quantification is called groundwater recharge estimation. In light of this challenge, studies and maps related to the provision of water resources are needed, with the purpose of indicating priority areas for conservation or restoration and to direct public policies to protect this natural resource. In Brazil, groundwater recharge studies that exist are predominantly focused on the river basin scale [8].

The incorrect management of river basins can generate serious threats to water availability, hindering surface and groundwater, because the dynamic of hydrologic systems is vulnerable to human actions [9,10]. An example is the reduced capacity for water infiltration into the soil observed in areas that are occupied with civil constructions or impermeable pavements [10,11]. Other problems, such as the intensification of erosion or floods, can occur in the development of soil permeability declines. In rural catchments, the inadequate uses of the land, among other factors, can amplify the risk of flood occurrence [12] reducing the levels of recharge.

Changes in the dynamics of hydrologic systems are particularly worrying in aquifer recharge zones, which are regions that enable water infiltration and percolation toward an aquifer system, which is defined as the geological system capable of storing and distributing a significant amount of water [13–15]. Moreover, recharge zones can be defined as areas where the soil surface favors the water infiltration and percolation [16,17]. Water can also be retained in the soil and slowly recharge the aquifer [13,18,19].

Some recharge zones are more efficient than others and for that reason are called preferential groundwater recharge zones [20]. The environmental protection of these special areas is important to conserve the quality and quantity of water resources. Thus, detailed information on the groundwater recharge process can aid better land use and cover distribution, and indicate the best areas for agricultural activities with the lowest groundwater contamination risks caused by the release of substances with high polluting potential, such as pesticides.

Despite the importance of sustainable management of aquifer recharge zones [7,21–23], in the state of Minas Gerais, Brazil, this topic has not yet been properly studied. Therefore, a better comprehension about the factors that affect groundwater recharge is necessary, as well as the mapping of recharge areas that consider the sustainable management potential of the basin. These evaluations should be robust and contain physical–environmental factors [10], such as soil characteristics, geology, vegetation cover, climate, and topography. When all these data are evaluated together, they allow for sustainable water use, meaning a use without compromising groundwater recharge [10]. Besides, under these circumstances the volume of water withdrawn from the aquifer system can be defined according to its natural capacity [24].

Numerous groundwater recharge estimation methods exist. However, they all entail some level of uncertainty [25]. In general, the practical and conceptual limitations of recharge estimation models occur because the available hydrological and hydrogeological data are sparse or fragmented, and because the spatial and temporal variations in recharge are significant [26]. This difficulty is significant for semi-arid areas [27], causing recharge estimations in these areas to be even more challenging. There are direct and indirect methods to evaluate the groundwater recharge potential. The direct methods include geological and geophysical explorations, gravimetrical and magnetic models, and perforation tests [10]. The indirect methods include hydrological and hydrogeological models [28,29], using geographical information systems (GIS) combined with field work [30,31]. Other studies have employed different methods to estimate groundwater recharge, which are comprised of tracer methods, water table fluctuation models, lysimeter methods and simple water balance techniques. Some of

these studies have used numerical groundwater models or dynamically linked them to hydrological models to estimate recharge variations under different climate and land cover conditions [32–36]. For example, Döll (2008) modeled global groundwater recharge using the WaterGAP Global Hydrological Model (WGHM), which has failed to reliably estimate recharge in semi-arid regions [37]. In that study, the influence of vegetation was not taken into account, even though many studies have showed the importance of this variable for estimating the groundwater recharge [32,38–42]. Moreover, Chowdhury et al. (2010) delineated groundwater recharge zones in West Medinipur district, India, using a GIS approach mixed with remote sensing and multi-criteria decision making techniques [22]. The input variables considered in that study were geomorphology, geology, drainage density, slope and aquifer transmissivity. In general, the choice of a method should consider the precision level needed, the project execution viability, and the available financial resources.

Among methods available in the literature, Costa et al. (2019) proposed one for the evaluation of groundwater recharge potential based on the water balance approach that considers climatic variables, water runoff, and the percolation of water into the soil profile [10]. The authors obtained results for mean annual recharge similar to those calculated by the hydrograph recession curve analysis, which has been used as a validation method. Furthermore, they identified areas with larger recharge potential and suggested management practices to improve groundwater recharge in those areas, making their study a valuable tool for the sustainable use of groundwater and protection of recharge areas [10].

Despite the positive results obtained by Costa et al. (2019), it is worth noting that the land management practices were a consequence of groundwater recharge assessments, not a contributing factor to groundwater recharge included in the model. Indeed, all parameters included in Costa's water balance model were physical, while land management practices had no role, regardless of their potential to dynamically affect groundwater recharge [10]. Thus, the coupling of physical factors and land management practices in a recharge estimation model could be a motivation (and a novelty) for a subsequent study. Before the publication of Costa et al. (2019), Costa et al. (2017) [43] conducted a study in Minas Gerais and developed a method based on multi criteria analysis, which was efficient to map a so-called Conservative Use Potential (PUC). The PUC method weights a considerably large number of variables considering their importance for sustainable land use, including several variables linked to land management practices, such as drainage, soil depth and fertility, erosion potential, and land capability. Hence, one possible route to realize our research motivation would encompass including the PUC, as determined by Costa et al. (2017) [43], within the framework of the Costa et al. (2019) groundwater recharge method [10].

The general purpose of this study is therefore to take that step forward and embed the concept of PUC in the groundwater recharge method of Costa et al. (2019) [10]. In that method, a parameter is defined to measure water percolation based on the soil's effective porosity and hydraulic conductivity. The method presented in this work replaces porosity by three parameters strongly influenced by management practices, which are soil texture, drainage, and profile depth. The replacement has the specific purpose to check whether this set of variables responds more effectively to land use changes than the original variable (porosity). The model was tested in the Pandeiros River Basin (PRB), Brazil, to generate a spatially explicit map of groundwater recharge potential, at regional scale (1:100.000). This map has the potential to subsidize indications of preferred areas for restoration, recovery, and protection, ensuring a more sustainable water resources management in this basin.

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

#### *2.1. Study Area*

The Pandeiros River basin (PRB) is located in the northern state of Minas Gerais, Brazil, and has an area of 396,028 ha, which encompasses part of Januária, Bonito de Minas, and Cônego Marinho municipalities (Figure 1). The climate of the region is predominantly dry with mean annual temperature of 24.6 ◦C, which occasionally reaches a maximum of 33 ◦C in October, and a minimum of 14 ◦C in July [44]. Rainfall is concentrated in April to September, since the climate is semiarid, with mean annual rainfall of approximately 1,050 mm year [45–47]. The longest and largest tributaries of the Pandeiros River are the Catolé, Suçuarana, Borrachudo and Macaúbas streams.

**Figure 1.** Location of the Pandeiros River basin in the state of Minas Gerais, Brazil.

According to the Brazilian Institute for Geography and Statistics (2018), the municipalities comprising the PRB have a total population of 86,311 inhabitants, most of them living in Januária. The largest part of this population live in rural areas, especially in the municipalities of Bonito de Minas and Cônego Marinho (Table 1) [48].


**Table 1.** Area and population living in the municipalities of Pandeiros River basin. Source: Brazilian Institute for Geography and Statistics. The population numbers refer to 2018.

The PRB is in an area with transitional vegetation, presenting phytophysionomies of Cerrado and Caatinga biomes [49]. This ecotone is characterized by swamp regions that contain the springs of the São Francisco River. These springs are responsible for the reproduction of most fishes that live between the Três Marias (MG) and Sobradinho (BA) dams [50].

The relief is predominantly flat. The plain was formed by the filling of the São Francisco Depression with sediments that were sourced from the erosion of rocks from the São Francisco Plateau [51]. This process is also responsible for a small proportion of steep slope areas in the basin, as shown in Figure 2 and quantified in Table 2 [52].

**Figure 2.** Slope map of the Pandeiros River basin.


**Table 2.** Area of slope classes in the Pandeiros River basin.

As regards the local geology, there is the presence of alluvial deposits over the main drains in the area and in wetland regions (Veredas), which are the results of natural and anthropogenic erosion processes with consequent transport and deposition of sediments. Figure 3 and Table 3 show the spatial distribution and proportion of lithotypes in the studied basin, respectively.

**Figure 3.** Lithological map of the Pandeiros River basin.


**Table 3.** Spatial distribution and proportion of lithotypes in the Pandeiros River basin.

The Soil Map of Minas Gerais [53] presents an area predominantly composed of red–yellow latosols, which cover more than 87% of the hydrographic basin, followed by fluvic neosols (5.48%). The occurrence of quartzarenic neosols, litolic neosols, cambisols and melanic gleysols are related to much smaller areas (Figure 4 and Table 4).

**Figure 4.** Soil map of the Pandeiros River basin.


**Table 4.** Spatial distribution and proportion of soil classes in the Pandeiros River basin.

The spatial distribution of the land use and cover classes in the basin (Figure 5) shows the predominance of a typical Cerrado vegetation (savanna) of low to medium size [54], which is usually associated with the occurrence of latosols and sandstone regions. These are found mainly in the central part, towards the northern and northeast of the Pandeiros River basin. This phytophysiognomy covers 183,719.88 ha in the basin, representing 46.3% of its area (Table 5). The second most significant land use and cover class in the study area is the dense Cerrado (21.5%), followed by the sparse Cerrado (10.42%).

Total 396,028.29 100

**Figure 5.** Land use and cover map of the study area, Pandeiros River basin.



#### *2.2. Material*

The materials used in this study (Table 6) consisted of (i) a digital elevation model (ALOS PALSAR), with a spatial resolution of 12.5 meters; (ii) a land use and cover map (scale of 1: 25.000); (iii) the soil map of Minas Gerais state (scale of 1: 600.000); (iv) values of groundwater recharge calculated by the PUC method and the hydraulic conductivity of each soil class in the basin; (v) rainfall and evapotranspiration data from meteorological stations located in municipalities near the basin; and vi) flow data records from the station code 45,250,000 located near the mouth of the Pandeiros River.



#### *2.3. Methodology*

The spatially explicit groundwater recharge potential in the Pandeiros River basin was evaluated using the methodology proposed by Costa et al. (2019) [10], with the adjustments described below. The workflow was divided into 5 main steps: (i) acquisition of a land use and cover map, digital elevation model, soil type map, and climate data (rainfall and evapotranspiration); (ii) calculation of surface runoff using the slope length and steepness factor, and runoff coefficients for the land use and cover types; (iii) calculation of water percolation based on the PUC method [43] (replacing the effective porosity of Costa's approach and representing the proposed methodological improvement) and adapted hydraulic conductivity values with fuzzy logic [55,56]; (iv) calculation of groundwater recharge in different points of the basin and a mean value for the whole basin, using a geographical information system; (v) validation of results by comparing the previously calculated mean groundwater recharge with the value estimated by the hydrograph recession curve analysis (Figure 6).

**Figure 6.** Flowchart for the groundwater recharge potential calculation and model validation. Adapted from Costa et al. (2019) [10]. Abbreviations: LS—slope length and steepness factor [57]; C—runoff coefficient (Table 7); RF—runoff factor (Equation (1)); PUC—Conservative Use Potential [43]; Ks—soils' hydraulic conductivity (Table 8); PF—water percolation factor (Equation (2)).

In the first step, the mean annual rainfall and evapotranspiration of municipalities near the study area were estimated using data from the Brazilian Institute for Meteorology (INMET) relative to the 2009–2018 period, which were obtained from meteorological stations in Arinos (MG), Januária (MG), Montes Claros (MG), Salinas (MG), Cariranha (BA), Espinosa (MG), Formoso (MG), Posses (GO), and Brasília (DF). Longer temporal series would be more adequate for estimating groundwater recharge. However, these records were not available. The flaws in the temporal series were resolved by using the regional weighting method, and the information was interpolated through inverse distance weighting (IDW) raised to the power of two [58]. This method was used because it raises the importance of closer stations in the interpolation. The data were spatialized and trimmed up to the study area limits.

In the second step, the land use and cover map and the topographical information from the digital elevation model (slope length and steepness factor—LS factor) were used to calculate the surface runoff factor, based on the method proposed by Böhner and Selige (2002) [59]. The runoff factor was calculated according to Equation (1), reproduced from Costa et al. 2019 [10].

$$RF = 1 - \left(\mathbb{C} + L S\_{\text{FULZZ}Y}\right) \tag{1}$$

where *RF* is the runoff factor (dimensionless); *C* is the runoff coefficient (dimensionless values adopted from Table 7); *LSFUZZY* is the slope length and steepness factor estimated using the method of Desmet and Govers (1996) [57], but recast to the 0–1 range using a fuzzy logic algorithm (the steeper the slope the closer to 1).

**Table 7.** Runoff coefficients for each land use and cover class in the Pandeiros River basin. Adapted from the ASCE—American Society of Civil Engineers, 1969, and Costa et al., 2019 [10,60].


In the third step, the water percolation was calculated considering the soil classes in the basin [53]. A systematic literature review was performed to set up hydraulic conductivity values for each soil class, at the second category level of Brazilian soil classification system (SiBCS) [61]. The hydraulic conductivity values were established according to Freire et al. (2003), Costa et al. (2019), Pedron (2011) and Amaral (2017) [10,55,56,62].

**Table 8.** Soil classes and respective *KsFuzzy* values. The scores of groundwater recharge parameter were set up by the PUC method in the Pandeiros River basin [10,55,56,62].


Different from the method proposed by Costa et al. (2019) [10], the water percolation factor in this model was calculated using the groundwater recharge category of the PUC method, developed by Costa et al. (2017) [43]. The PUC is a method that allows for mapping areas of a basin based on their limitations and potentialities for conservationist land use, through the combined assessment and weighting of several environmental variables (soils, geology, and geomorphology) [43].

The PUC assigns values from 1 to 5 to the different classes of lithology, slope and soil in the watersheds from the state of Minas Gerais, Brazil. The analyses are focused on groundwater recharge, agricultural use potential and soil resistance to erosion in the catchments [43]. For the identification of lithologies, slopes and soil classes existing in Minas Gerais, the available official databases are used.

The attribution of grades for the different types of lithologies took into consideration their potential to provide nutrients (greater weight for rocks with higher absolute content of essential macro elements to plants) and their susceptibility to weathering processes (considering the main mineral constituents and scored according to their resistance to weathering, based on Goldish's stability [63]). Regarding the slope parameter, the same weights were attributed for groundwater recharge potential for agricultural use and resistance to erosion. For groundwater recharge, it was considered that the slope has a direct relationship with the flow velocity and the opportunity time for water infiltration. The higher the slope, the higher the water velocity and the shorter the time for water infiltration. In this context, the mountainous relief received a weight 1 and the flat relief a weight 5.

For the attribution of grades to different soils, the variables texture, drainage, effective depth and fertility were considered. For groundwater recharge, the attribute fertility was disregarded and the classes of soils characterized as favorable to water infiltration and percolation received greater weight. The recharging potential for each soil type was obtained by the simple average of the values of texture, drainage and effective depth, normalized so that the final scale was in the range of 1 to 5.

In this work, only the soil parameter of groundwater recharge was used, which assigns the basin´s soil classes a score from 1 to 5. This parameter takes into account the effective depth, texture and drainage of each soil class at the first category level of SiBCS [61]. Thus, to the soil that presents considerable effective depth, satisfactory drainage and texture favorable to infiltration, a higher score regarding groundwater recharge is attributed.

However, in order to adapt the PUC method to the recharge model, a rescaling process was implemented, whereby the ratings from 1 to 5 were recast to the range 0 to 1. The values used in the model are represented in Table 8. The percolation factor was evaluated using Equation (2), reproduced from Costa et al. 2019 [10]:

$$PF = PUC \times Ks\_{FULZY} \tag{2}$$

where *PF* is the water percolation factor (dimensionless); *KsFUZZY* (dimensionless) is the soil's hydraulic conductivity fitted to the 0 to 1 range by the fuzzy logic algorithm (the higher the soil's hydraulic conductivity in the class the closer to 1); *PUC* are the scores of groundwater recharge set up by the PUC method but fitted to the 0 to 1 range.

In the fourth step, the groundwater recharge potential was calculated for each point in the basin, using Equation (3), according to Costa et al. 2019 [10]:

$$RPot = \left[ (P - ET\_R) \times RF \times PF \right] \times 10 \tag{3}$$

where *RPot* is the groundwater recharge potential (m<sup>3</sup> ha−<sup>1</sup> year<sup>−</sup>1); *P* is the mean annual rainfall depth (mm year<sup>−</sup>1); *ETR* is the mean evapotranspiration (mm year<sup>−</sup>1); *RF* is the surface runoff factor; and *PF* is the water percolation factor.

The results were validated in the fifth step through comparison of the calculated mean groundwater recharge with a homologous median value estimated by the hydrograph recession method based on the Maillet equation [64,65]. The hydrograph for the Pandeiros River was drawn from daily average stream flow data measured at the hydrometric station nº 4425000 located in the Pandeiros River mouth (Pandeiros River Dam), compiled from the Hidroweb portal [66]. Streamflow data from 2013 to 2018 were the most recent and more continuous in the historical series, showing no flaws or absent data. Thus, these data were used to calculate the Maillet equation (Equation (4)). To separate and analyze the recession curve and the recession days, the methodology proposed by Barnes (1939), Dewandel et al. (2003) and Kovacs et al. (2005) was used [64,65,67–69]:

$$Q\_T = Q\_0 \times e^{-\alpha t} \tag{4}$$

where *QT* is the flow at time *t* (m3 s<sup>−</sup>1); *Q0* is the flow at the beginning of a recession (m<sup>3</sup> s<sup>−</sup>1); α is the coefficient of recession; *t* is the time (days); and *e* is the basis of Neperian logarithm (2.71828).

Thus, the coefficient of recession can be determined numerically, based on the logarithmic form of Equation (4), represented and rearranged in Equation (5):

$$\alpha = \frac{\text{LogQ}\_0 - \text{LogQ}\_t}{0.4343t} \tag{5}$$

Subsequently, the groundwater recharge volume was calculated using Equation (6):

$$V = \frac{Q\_0 \times t'}{a} \tag{6}$$

where *V* is the recharge volume (m3); *Q0* is the flow at the beginning of recession (m<sup>3</sup> s−1); *t'* is the converter of the *t* unit (days into seconds; 86,400); α is the coefficient of recession (dimensionless).

The constant of recession (α; Equation (5)) is dependent on the aquifer characteristics and therefore should not vary significantly from year to year. On a *Q* versus *t* plot (hydrograph), where the *Q* values are represented in logarithmic scale and the values are in linear scale, the baseflow within a hydrologic year (from the recharge period to the end of recession) should define as a straight line, the slope to which is related in terms of α. If *t*cycle is the time of a log cycle for discharge, meaning the time for discharge to change from 1 to 10 m3/s, from 10 to 100 m3/s, and so forth, then *LogQ*<sup>0</sup> <sup>−</sup> *LogQt* = 1 and α = 1/(0.4343*t*cycle). This simplified representation of Equation (5) is frequently used in the calculation of α and will be adopted in the present study. Conversely, the values of *Q*<sup>0</sup> can vary in response to the annual variations of precipitation. In this case, a value of *Q*<sup>0</sup> should be calculated for each hydrologic year, while mean ± standard deviation values are derived therefrom.

#### **3. Results**

The interpolation of precipitation data from 2009 to 2018, obtained from climatologic stations close to the Pandeiros River basin, showed mean rainfall depths ranging from 904.7 to 1056.3 mm year−<sup>1</sup> (Figure 7).

**Figure 7.** Spatial distribution of rainfall in the Pandeiros River basin.

The evapotranspiration ranged from 729.6 to 856.5 mm year−<sup>1</sup> (Figure 8). The spatial distribution of this variable was similar to rainfall; the lowest values were found in the southeast, east, and northeast regions, whereas the highest values were found in the northwest region of Pandeiros River basins.

**Figure 8.** Spatial distribution of evapotranspiration in the Pandeiros River basin.

The evaluation of the runoff factor exposes a strong effect of runoff coefficients on the RF, as observed when the land use and cover map (Figure 5) is compared with the runoff map (Figure 9). The map of Figure 9 showed a lower groundwater recharge in urban areas, which have a higher runoff. In contrast, higher potential for groundwater recharge takes place in dense vegetation areas, indicating that the presence of vegetation decreases surface runoff.

**Figure 9.** Spatial distribution of runoff factor in the Pandeiros River basin.

The water percolation factor ranged from 0.001 to 0.7 (dimensionless), expressing the combined variation of hydraulic conductivity values and PUC scores fitted to the 0–1 range (Figure 10 and Table 8). Areas with melanic gleysols and fluvic neosols presented the lowest water percolation factors, decreasing the groundwater recharge in these areas because of their low hydraulic conductivity.

**Figure 10.** Spatial distribution of water percolation in the Pandeiros River basin.

The map of Figure 11 shows the groundwater recharge potential of the Pandeiros River basin ranging from 0 to 122.7 mm year−1, with a mean value of 93.99 mm year−1. The areas with higher groundwater recharge potential are located in regions with dense tree vegetation cover, areas with flat or slightly wavy relief, and areas with developed and structured soils, where the porosity and hydraulic conductivity allow water percolation to the water table. These areas are distributed throughout the basin and are found in the three municipalities that encompass the basin.

The areas with lowest groundwater recharge potential are represented in Figure 11 with a red color and are located in the urban area and in areas with the presence of melanic gleysols and fluvic neosols. The reduced potential is due to the soil sealing process occurring in urban areas, and to low hydraulic conductivities that decrease water percolation and consequently recharge in the aforementioned soil types.

The validation of PUC-based recharge estimates is depicted in Table 9. The mean groundwater recharge is 93.99 mm year−<sup>1</sup> which is close to the median value obtained with the hydrograph recession method (87.2 mm year−1). The difference between the two values is just 6.6 mm year−<sup>1</sup> or 7.3%. The hydrograph used to calculate groundwater recharge based on base flows is shown in Figure 12. Hydrographs describe the succession of peaks representing the watershed response to a precipitation event, which are separated by baseflow segments that describe the aquifer response to drainage. Hydrograph recession curves are usually separated in the quick flow stage, which depicts the runoff and water infiltration towards the saturated zone, and the baseflow stage when only the saturated zone discharges [68]. The baseflow discharge is the most representative feature of an aquifer's global response because it is less influenced by the temporal and spatial variations on infiltration [69]. Generally, hydrographs are analyzed together with rainfall data. The peak of a hydrograph rising curve shows the highest values of stream flows, which take place in the months when rainfall values

are the highest. In these periods the superficial runoff also reaches its highest values and it decreases during the recession time, marked by the red segments in Figure 12.

ϭ ϭϬ ϭϬϬ ϭϬϬϬ ϬϭͲϬϭͲϮϬϭϯ ϬϭͲϬϯͲϮϬϭϯ ϬϭͲϬϱͲϮϬϭϯ ϬϭͲϬϳͲϮϬϭϯ ϬϭͲϬϵͲϮϬϭϯ ϬϭͲϭϭͲϮϬϭϯ ϬϭͲϬϭͲϮϬϭϰ ϬϭͲϬϯͲϮϬϭϰ ϬϭͲϬϱͲϮϬϭϰ ϬϭͲϬϳͲϮϬϭϰ ϬϭͲϬϵͲϮϬϭϰ ϬϭͲϭϭͲϮϬϭϰ ϬϭͲϬϭͲϮϬϭϱ ϬϭͲϬϯͲϮϬϭϱ ϬϭͲϬϱͲϮϬϭϱ ϬϭͲϬϳͲϮϬϭϱ ϬϭͲϬϵͲϮϬϭϱ ϬϭͲϭϭͲϮϬϭϱ ϬϭͲϬϭͲϮϬϭϲ ϬϭͲϬϯͲϮϬϭϲ ϬϭͲϬϱͲϮϬϭϲ ϬϭͲϬϳͲϮϬϭϲ ϬϭͲϬϵͲϮϬϭϲ ϬϭͲϭϭͲϮϬϭϲ ϬϭͲϬϭͲϮϬϭϳ ϬϭͲϬϯͲϮϬϭϳ ϬϭͲϬϱͲϮϬϭϳ ϬϭͲϬϳͲϮϬϭϳ ϬϭͲϬϵͲϮϬϭϳ ϬϭͲϭϭͲϮϬϭϳ ϬϭͲϬϭͲϮϬϭϴ ϬϭͲϬϯͲϮϬϭϴ ϬϭͲϬϱͲϮϬϭϴ ϬϭͲϬϳͲϮϬϭϴ ϬϭͲϬϵͲϮϬϭϴ ϬϭͲϭϭͲϮϬϭϴ 6WUHDPIORZPñV 7LPH'D\V**Ϳ** YϬ Y YϬ YϬ <sup>Ϭ</sup> YϬ YϬ ƚĐLJĐůĞсϳϴϵĚĂLJƐ

**Figure 11.** Spatial distribution of groundwater recharge potential in the Pandeiros River basin.

**Figure 12.** Hydrograph of Pandeiros River basin, state of Minas Gerais, Brazil, used to calculate the recession constant value.


**Table 9.** Groundwater recharge in the Pandeiros River basin, estimated by the spatially distributed PUC-based and hydrograph recession analysis methods.

In the assessment of recharge using the method of hydrograph recession analysis, we looked for straight line segments corresponding to base lows. These segments allowed us to draw the corresponding fitting lines (dashed red lines). As would be expected, the lines are all virtually parallel because the line slope is solely dependent on the aquifer characteristics and dimension, which are invariant at the timescale of a few years (present case). For all years, the point in the graph where the fitting line intercepted the hydrograph at the upper flows (left edge point) was defined as the Q0,t0 point, i.e., the point where the recession period began. On average Q0 = 10.8 <sup>±</sup> 1.8 m3 s<sup>−</sup>1. The coefficient of recession was estimated by the simplified version of Equation (5), α = 1/(0.4343*t*cycle). In the Pandeiros River basin, the estimated *t*cycle is 789 days, and therefore α = 0,002918. Using Q0 and α in Equation (6) results in the recharge value of *<sup>V</sup>* <sup>=</sup>98.8 <sup>±</sup> 20.3 mm year<sup>−</sup>1.

Figure 13 shows the spatialization of groundwater recharge within the municipalities that constitute the studied basin. Although the municipality of Januária encompasses a large area of the Pandeiros River basin, approximately 212,956.8 ha, this municipality presents a higher proportion of melanic gleysols, fluvic neosols, and litolic neosols, which are soils with the lowest water percolation factor (PF) assessed by the model. Moreover, it presents a considerable proportion of exposed soil, sparse Cerrado vegetation, and cultivated soils in the northwest region of the basin, which increases the runoff potential. For these reasons, this municipality presented a mean annual groundwater recharge of 90.30 mm year<sup>−</sup>1.

The municipality of Cônego Marinho presented the highest mean annual groundwater recharge (105 mm year<sup>−</sup>1). This municipality occupies a smaller area of the basin (approximately 26,083.86 ha) than the other municipalities, but it is in a region where soil and topography favor groundwater recharge.

Bonito de Minas presented mean annual groundwater recharge of 97.18 mm year<sup>−</sup>1. The area of the basin within this municipality is approximately 156,987.95 ha. Despite the considerable presence of fluvic neosols and melanic gleysols in this region of the basin, the proportion is lower than that of Januária. The urban area of Bonito de Minas within the basin is small, representing 0.03% of the total area of the basin, which does not significantly affect the mean annual groundwater recharge.

The land use and cover map (Figure 5) and the groundwater recharge map (Figure 11) showed that the areas with forest cover presented the best recharge values. This type of land use and cover combined with flat or slightly wavy relief and areas overlaid with red–yellow latosols proved preferable for groundwater recharge.

Contrastingly, regions presenting such land use and cover as urban areas, exposed soil, ground vegetation, poorly structured soils and low hydraulic conductivity (fluvic neosols, litolic neosols, haplic cambisols, melanic gleysols), combined with steeper areas, presented lower recharges.

**Figure 13.** Spatial distribution of groundwater recharge potential in the municipalities of Pandeiros River basin.

#### **4. Discussion**

The rainfall depths (Figure 7) are consistent with Neves (2011), who evaluated a historical series (1931–1990) in the same region and reported a mean rainfall depth of 1050 mm year−<sup>1</sup> [45]. However, the Brazilian Institute for Geography and Statistics (2014) considers values below 800 mm year−<sup>1</sup> to define the semiarid region of Brazil [70], denoting that the rainfall depths of the study area were relatively high for a semiarid region.

Yang et al. (2016) showed that the rainfall depth can be a determinant for evapotranspiration and temperature values [71]. Therefore, practically all water that leaves the system (output) in regions with dry climate and low rainfall, such as the Brazilian semiarid regions, is due to evapotranspiration [72]. According to Loos, Gayler, and Priesack (2007), evapotranspiration is one of the most critical variables, with a high impact on water loss in similar regions [73].

The spatial distribution of runoff factor showed that cover plants are essential to maintaining the water cycle and to protect the soil against the impacts of raindrops. Moreover, the presence of vegetation increases soil porosity and permeability by the action of roots, thereby reducing runoff, and keeping soil moisture in the vicinity of organic colloids [74].

Soils under forest are characterized by expressive plant residue layers (litter fall) and by an A horizon rich in organic matter, which enables a higher soil aggregation, preserving its porosity [75]. Soils under forest usually present significant porosity, mainly macropores due to dead roots and animal holes, which are important to facilitate water infiltration and recharge. Therefore, water infiltration capacity is usually more expressive in areas with forest vegetation [75,76] than in pastureland or cropland, as found in the present work, resulting in a lower surface runoff.

Costa et al. (2019) [10] used a groundwater recharge model and reported a higher surface runoff in urban areas than in vegetated areas. Urban areas have drainage systems, street pavements, and infrastructures that hinder infiltration of rainwater into the soil, decreasing the recharge potential [10].

Comparing the first model employed by Costa et al. (2019) with the present application, it is interesting to see how including the groundwater recharge parameter from the PUC method improved the spatialization of groundwater recharge, mainly because the first model [10] used solely total porosity values for each soil class in the Jequitiba River basin [10]. This parameter is not the best to assess drainage and water percolation, because according to Silva et al. (2013), the micro and macro porosity better reflect the movement of water through the soil profile [77]. Moreover, the PUC method uses a management approach that was not taken into consideration in the PF in the first model. Another important difference from the first model to the present one is the separation until the second category level of the neosols class, which has shown more reliable spatialization and values for the different category levels of this soil class. The differences between the physical characteristics (structure, texture and effective depth) of fluvic neosols, litolic neosols and quartzarenic neosols) make it questionable to give the same values of hydraulic conductivity and porosity to this soil class.

Regarding the validation stage of the work, the higher the slope of the recession curve, the faster the depletion of the water table reserves and the higher the demand of this system to regulate the surface system; whereas the lower the slope, the lower the time of depletion of water resources. The surface and groundwater form a system with mutual contribution, thus, any change in one will affect the other in the short or long term.

The comparison of the methods (spatialization; recession curve analysis) showed a small difference (5.2%) because of their particularities, since each method has inherent and inevitably uncertainty levels [78]. The spatialization method probably has a higher uncertainty in groundwater recharge values because the number of parameters involved with the calculations is much larger [10]. However, the recession curve analysis provides only an average groundwater recharge estimate (mean) for the entire basin, whereas the spatialization method provides one estimate for each point in the basin.

The small difference between the results obtained with the spatially distributed PUC-based and hydrograph recession methods attributes to geology a limited role in the estimation of recharge in the Pandeiros River basin, and that this process is predominantly controlled by the infiltration capacity and profile depth of the soil. The effect of vadose zone thickness on the recharge was recently reported in China [79]. Information on the characteristics of a river basin and its potentialities and limitations are essential to an adequate management of water resources [10,80,81].

The potentialities found for the municipality of Cônego Marinho include the prevalence of red–yellow latosols, which are soils that favor groundwater recharge because of their structural characteristics, such as the occurrence of macropores that increase hydraulic conductivity [43]. Moreover, the region has few cultivated areas and has a predominance of typical Cerrado vegetation, which also favor groundwater recharge, according to the model. Thus, the mean annual groundwater recharge was higher in this municipality than in the other two in the basin.

The considerable overlay of fluvic neosols and melanic gleysols in Bonito de Minas is compensated for by a high proportion of typical Cerrado and dense Cerrado vegetation and by the presence of few cultivated areas, which decreased the runoff potential: Consequently, the mean annual water recharge in this municipality was higher than that in Januária. The presence of native vegetation favors the groundwater recharge, because the forest reduces the surface runoff, favors water percolation, and maintains the soil physical and mechanical stability, assisting in the storage of water and the supply of groundwater [10]. Several studies evaluated the effect of land use and cover and the benefits of forested areas to groundwater recharge [42,82–88].

Considering these results, the method used in the present work provides aid in the better management of river basins by considering their potentialities and limitations, not generalizing a mean value for the whole area evaluated. Therefore, preferred areas can be identified for recharge, thus directing public policies and conservationist actions for each area, according to their needs.

#### **5. Conclusions**

The groundwater recharge potential was higher in areas covered with forests, located in plains or slightly wavy relief areas, or overlaid with red–yellow latosol. These are the areas to protect in the watershed management plan if recharge is to be favored or restored. The soil classes and their structural attributes, as well as the land use and cover types were considered as key factors for groundwater recharge. The results showed that areas with higher groundwater recharge potential were concentrated in the municipality of Cônego Marinho, followed by Bonito de Minas, and Januária. Areas with a presence of melanic gleysols and fluvic neosols presented the worst responses in the model.

As made evident from the results, this study should be used as a tool for the management of water resources in the Pandeiros River basin, because the preferred recharge areas could be successfully identified. It is urgent therefore that public policies and conservationist actions are enforced in these areas to improve natural groundwater recharge, and hence increase the accessible water volume to the local population. The adjustment of irrigation methods, adoption of soil preservation practices to improve water infiltration, seasonal storage of surface water in areas of low recharge potential and the preservation of forest vegetation, are examples of feasible actions. Moreover, this work provided subsidies for further studies that seek methods for the spatialization of groundwater recharge potential in river basins with a key role assigned to management practices.

**Author Contributions:** Conceptualization, M.A.T., M.S.d.M. A.M.d.C., P.K.M., J.H.M.V.; methodology, M.A.T., M.S.d.M., A.M.d.C. and J.H.M.V.; validation, M.A.T and F.A.L.P.; resources, A.M.d.C.; writing—original draft preparation, M.A.T. and M.S.d.M.; writing—review and editing, M.A.T., M.S.d.M., A.M.d.C., P.K.M., F.A.L.P and L.F.S.F.; supervision, A.M.d.C.; funding acquisition, A.M.d.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the research project "Sustentabilidade da Bacia do Rio Pandeiros" (Sustainability of Pandeiros River Basin), sponsored by grant APQ-03773-14 from FAPEMIG – Fundação de Amparo à Pesquisa do Estado de Minas Gerais that included a research scholarship for the author Maíse Soares de Moura. The manuscript translation from Portuguese to English was financed by the CAPES – Coordenação de Aperfeiçoamento de Pessoal de Nível Superior. For the author integrated in the CITAB research centre, the research was financed by National Funds of FCT–Portuguese Foundation for Science and Technology, under the project UIDB/04033/2020. For the author integrated in the CQVR, the research was supported by National Funds of FCT–Portuguese Foundation for Science and Technology, under the projects UIDB/00616/2020 and UIDP/00616/2020.

**Acknowledgments:** This study was conducted within the work plan of research group "GEISS – Grupo de Estudos Integrado em Solos e Sustentabilidade", in its research line on "Gestão Sustentável de Recursos Hídricos e Segurança Hídrica". We wish to thank the support of the team of Soil and Environment Laboratory of the Federal University of Minas Gerais for the help in the analyses made in the study.

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

#### **Nomenclature**

The list of mathematical symbols and their measurement units as used in the present article are listed below in alphabetical order:


#### **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/).

## **The Assessment of Hydrological Availability and the Payment for Ecosystem Services: A Pilot Study in a Brazilian Headwater Catchment**

**Mariana Bárbara Lopes Simedo 1,6, Teresa Cristina Tarlé Pissarra 1,6, Antonio Lucio Mello Martins 2, Maria Conceição Lopes 2,6, Renata Cristina Araújo Costa 6, Marcelo Zanata 3,6, Fernando António Leal Pacheco 4,6,\* and Luís Filipe Sanches Fernandes 5,6**


Received: 24 August 2020; Accepted: 27 September 2020; Published: 29 September 2020

**Abstract:** The assessment of water availability in river basins is at the top of the water security agenda. Historically, the assessment of stream flow discharge in Brazilian watersheds was relevant for dam dimensioning, flood control projects and irrigation systems. Nowadays, it plays an important role in the creation of sustainable management plans at the catchment scale aimed to help in establishing legal policies on water resources management and water security laws, namely, those related to the payment for environmental services related to clean water production. Headwater catchments are preferential targets of these policies and laws for their water quality. The general objective of this study was to evaluate water availability in first-order sub-basins of a Brazilian headwater catchment. The specific objectives were: (1) to assess the stream flow discharge of first-order headwater sub-basins and rank them accordingly; (2) to analyze the feasibility of payment for environmental services related to water production in these sub-basins. The discharge flow measurements were conducted during five years (2012 to 2016), in headwaters in a watershed on the São Domingos River at the Turvo/Grande Watershed, represented as the 4th-largest hydrographic unit for water resources management—UGRHI-15 in São Paulo State, Brazil. A doppler velocity technology was used to remotely measure open-channel flow and to collect the data. The discharge values were obtained on periodic measurements, at the beginning of each month. The results were subject to descriptive statistics that analyzed the temporal and spatial data related to sub-basins morphometric characteristics. The discharge flows showed space–time variations in magnitude between studied headwater sub-basins on water availability, assessed based on average net discharges. The set of ecological processes supported by forests are fundamental in controlling and recharging aquifers and preserving the volume of water in headwater in each sub-basin. The upstream inflows influence downstream sub-basins. To avoid scarcity, the headwater rivers located in the upstream sub-basins must not consider basin area as a single and homogeneous unit, because that may be the source of water conflicts. Understanding this relationship in response to conservationist practices installed

uphill influenced by anthropic actions is crucial for water security assessment. The headwaters should be considered a great potential for ecosystem services, with respect to the "provider-receiver" principle, in the context of payments for environmental services (PES).

**Keywords:** flow; water discharge ecosystem services; payments for environmental services; land use; riparian forest

#### **1. Introduction**

The study of water availability in watersheds is fundamental for the demonstration of water potential and hydrological behavior of a region. It also helps with increasing the capacity of a population to safeguard access to adequate quantity and acceptable quality of water to sustain well-being and the environment. Surface and groundwater reserves are considered strategic and essential components of ecosystems and are fundamental for economic, social, and sustainable development [1,2]. The availability, quality, management, and governance of water resources are currently at the center of technical and scientific discussions [3–6].

Although most of the planet Earth's surface is occupied by water, 97.5% of available water is salty, and only 2.5% is fresh water. From the percentage of fresh water, 68.9% is concentrated in glaciers, polar ice caps, or mountainous regions, 29.9% in groundwater, 0.9% in other reservoirs, and only 0.3% make up the portion of surface fresh water present in rivers and lakes [7]. Brazil has 13% of the world volume of fresh water available, with a higher concentration in the Amazon region, where there is less population and less demand, according to the National Water Agency (ANA) [8]. The State of São Paulo, which is the most populous, has 1.6% of Brazilian fresh water [9]. The multiple uses of water resources are diversified—public supply, food production, hydroelectricity generation, navigation and industrial development, and their intensity is related to social, agricultural, and industrial development [1].

In the agriculture sector, approximately 70% of total fresh water is used in production activities. This consumption tends to increase with population growth and demand for food [10,11]. Thus, the need to obtain information and metrics that contribute to the establishment of planning measures and future management policies for water security is evident. Water security involves ensuring water quantity and quality acceptable for livelihood. It is related to the development of solutions to manage and mitigate impacts of scarcity and possible risks regarding environmental conditions and climate uncertainties, in order to guarantee human well-being, socioeconomic development and preservation of ecosystems [12–17].

The management of water resources in an integrated planning policy is increasingly needed to ensure water security in the future [3,15,16,18]. Studies showed that more than half of the global population has experienced severe water shortages for at least one month in a year, and climate change is affecting the reliability of supplies and infrastructures available in many regions [15,16,19]. According to scenario-based studies from the Fifth Assessment of the Intergovernmental Panel on Climate Change (IPCC), and data modeled on WaterGAP3, it is estimated that by 2050 urban water demand will increase by 80% [20].

Brazil experiences annual imbalances between water supply and demand. A remarkable episode was in 2013–2014, in the southeastern region, in which drought affected approximately 80 million people in Minas Gerais, São Paulo, and Rio de Janeiro States [18,21–23]. There were considerable problems for public water supply, food production, energy and navigation, causing the loss of 5000 jobs and millions of tons of non-transported materials [18]. Water crises generate significant socioeconomic and environmental impacts. Thus, the development of strategies for adapting to change that represent improvements and efficiency in future availability of water in hydrographic basins is welcome.

In water security research, the working territorial unit is a watershed. That is the unit where the land use and land cover are best managed [24–28]. Research studies on water availability [29–35] and flow distribution [36–42] within this unit contribute to manage water scarcity or stress, understand availability and demand, establish grants, identify environmental impacts, and even implement public policies [43–48].

The information on flow (discharge) metrics of Brazilian drainage networks is relevant for preparing policies, projects for new dams and flood control, as well as creating sustainable management plans within the scope of land use and land cover management in the hydrographic basin. The basin is a place where several interactions between the elements of an ecosystem occur, and thus it can be considered that this space has an ecosystem function considering water and nutrient cycle, in regulation of gases, in transfer of energy, and in the water infiltration/discharge relationship. The ecosystem services generated by these functions trigger a series of benefits, directly or indirectly, that humans can be appropriate.

A single ecosystem service can be a product of two or more functions, and a single function can generate more than one ecosystem service. Discharge flow controlling management practices can enhance environmental services or minimize the impact of human activities in a territory, or even promote economic incentives aimed at conserving ecosystems and increasing these services, which is of great value [49–51]. Ecosystem services promoted by maintenance and conservation of water resources have extremely high economic and social value. Valuation of those services and their diversity have gained scientific and economic attention, and environmental service benefits and control functioning of ecosystems must provide human well-being [52–55].

The Payment for Environmental Services (PES) instrument is considered an efficient program, as it rewards those who produce or maintain environmental services and encourages those who would not promote these services in the absence of monetary stimulus [53,55,56]. The monetary compensation and reduction of tax charges can be applied, or even, use of public resources from municipal and state funds, charging water users, environmental compensation, and carbon credit [54,57,58].

In Brazil, PES initiatives are centered on projects related to water and watershed, carbon storage programs, biodiversity, and landscape protection. Some pioneering and developing examples in Brazil are: Water Producer Program from National Water Agency–ANA (https://www.ana.gov.br/ programas-e-projetos/programa-produtor-de-agua); Atlantic Forest Connection Project, financed by different development agencies (https://conexaomataatlantica.mctic.gov.br/cma/o-projeto/o-que-e); Conservative Water Project of Extrema, State of Minas Gerais (https://www.extrema.mg.gov.br); Ecocredit Program in the Montes Claros State of Minas Gerais (https://portal.montesclaros.mg.gov.br).

The assessment of hydrological availability in the Brazilian watershed on headwater sub-basins should therefore be recognized and promoted as a primary strategy to support the continued provision of ecosystem services in watersheds. The potential for ecosystem services provision, considering the principle of "provider-receiver", to establish PES must be investigated in the sequel. The costs to protect and manage those areas are substantial, and PES programs are a promising strategy. The landowner will receive a financial support as compensation to guarantee the provision of water ecosystem service to a willing buyer or beneficiary (e.g., Water Supplies Companies).

The general objective of this study was to evaluate water availability on sub-basins of first order and their contribution to water security assessment on Brazilian watersheds. The specific objectives were: (1) to assess the hydrological discharge of headwater sub-basins; (2) comprehend which first-order sub-basins are higher on water flow discharge; and (3) analyze the feasibility of payment for environmental services for the water-producing sub-basins.

#### **2. Study Area**

The experimental study area was the Olaria Stream Watershed, a tributary of the São Domingos river. This river contributes with flow to the Turvo/Grande watershed, the 4th-largest water resources management unit—UGRHI-15 in the State of São Paulo, Brazil (Figure 1), according to the Hydrographic Basins Committee of the Turvo and Grande Rivers—CBHTG (http://www.comitetg.sp.gov.br). The Olaria Stream Watershed covers 9.45 km<sup>2</sup> and is located between latitudes 21◦05 47 S and

21◦19 35 S and longitudes 49◦03 02 W and 48◦42 52 W Gr., expressed in UTM—Universal Transverse Mercator Projection System, Zone 22K, with altitudes varying from 507 to 616 m.

**Figure 1.** Olaria Stream Watershed, São Domingos River watershed, Turvo/Grande watershed, UGRHI-15, São Paulo State, Brazil.

The river flow supplies water to the municipalities of Catanduva, Catiguá, Tabapuã, and Uchoa, with a total of 169,232 inhabitants, according to data from the Brazilian Institute of Geography and Statistics—IBGE (https://cidades.ibge.gov.br). The area is socioeconomic-important for regional development, namely, the agricultural and industrial productive system. The sub-basin areas belong to Polo Regional Centro Norte, São Paulo Agribusiness Technology Agency—APTA (https: //www.apta.sp.gov.br), a Department of the State of Agriculture and Supply Secretariat, São Paulo, Brazil. This unit carries out agricultural research and experimentation with annual and perennial crops, in order to meet technological demand on agribusiness production chains.

The climate is classified as Cwa, defined as a subtropical dry winter (with temperatures below 18 ◦C) and hot summer (with temperatures above 22 ◦C), with an average annual air temperature of 20–22 ◦C [59]. The local water balance (Figure 2), according to 30 year data (1961 to 1990) from the Brazilian Agricultural Research Corporation—EMBRAPA [60], is characterized by: average temperature of 22.8 ◦C; total annual precipitation of 1388 mm; average annual precipitation of 116 mm; potential evapotranspiration of 1134 mm; soil water storage of 788 mm; real evapotranspiration of 1054 mm; water deficiency of 80 mm; and water surplus of 334 mm.

**Figure 2.** Water balance in the APTA-Pindorama, State of São Paulo, Brazil. Adapted from Brazil's climate database, available at the Brazilian Agricultural Research Corporation [60].

The watershed soil was classified as Red-Yellow Argisols, sandy/medium texture, with wavy and smooth wavy relief [61]. The predominant natural vegetation land cover was classified as semi-deciduous seasonal tropical forest, Atlantic Forest biome [62]. Land use was mainly distributed by urban areas, industries, agro-industries and agriculture, with a predominance of sugarcane and crop production systems such as citrus, rubber, grains, livestock, pasture, and eucalyptus [63]. The sub-basins uphill have an intense agricultural crop production and areas of native forests, considered as a Permanent Preservation Area (PPA) by Environmental Brazilian Laws.

#### **3. Materials and Methods**

Flow measurements in natural watercourses were carried out in order to determine the value of the surface runoff of a basin, its temporal variability, and the characteristics of the runoff. The activities of monitoring the amount of water in the Córrego da Olaria watershed in Pindorama-SP were initiated to understand the water in the hydrological processes and water production capacity of the sub-basins, obtaining stream flow data and the elaboration of quantitative analysis of flows, which will be indispensable for future actions and policies on land and water uses.

#### *3.1. Discharge Flow Monitoring in Sub-Basins*

The volume of water resources (e.g., discharge, flow) monitored in the Olaria stream watershed began in January 2012, as an activity of a research project entitled "Recovery of headwaters of the Polo Regional Centro Norte (Pindorama Experimental Station)". Subsequently, in 2013, it continued with the project "Monitoring water resources to assess changes associated to land use and soil management at Olaria stream watershed". The data acquisition was carried out in three sub-basins located at the headwaters of Olaria catchment and near its mouth where the stream debouches into the São Domingos river (Figures 3 and 4).

**Figure 3.** Drainage network of Olaria stream watershed, with identification of monitored sub-basins (1.1; 1.2; 2; 3).

**Figure 4.** Location of discharge flow monitoring points: (**a**) sub-basin 1.1; (**b**) sub-basin 1.2; (**c**) sub-basin 2; (**d**) sub-basin 3.

The uphill of sub-basin 1.2 (Figure 3) comprises the recovery and reforestation of a degraded area. This specific area had an aggravated erosion process due to inadequate soil management (Argisol—susceptible to erosion) (Figure 5). For decades, the coffee production system was installed uphill and later, pasture and cattle raising, with no soil conservationist practices [64]. The lack of vegetation for soil protection and its incorrect land use have severely degraded the sub-basin hills, resulting in a severe erosion process of a 700 m long gullet and some stretches, up to 15 m deep (Figure 5a) [64,65].

In order to contain the erosion process, recover the source and biodiversity, four uneven reservoirs were built (Figure 5b), interconnected by drainage channels (spillways on concrete stairs) (Figure 5c), which allowed water to pass through with a controlled flow. In addition, surrounding agricultural experimentation areas were managed with conservationist and maintenance treatments within natural vegetation land cover [64–66]. In 2011, reforestation of the surrounding stream was made within an Agroforestry System (AS), under a different management on the uphill of the sub-basin (Figure 5d). The purpose was to restore permanent preservation areas (PPA) and complement the other environmental recovery practices [66,67].

**Figure 5.** History of recovery and reforestation of degraded area in sub-basin 1.2: (**a**) intense erosion and gullies; (**b**) formation of four uneven reservoirs; (**c**) drainage channels; and (**d**) implementation of reforestation and the Agro-forestry System.

The selected points 1.1 and 1.2 belong to sub-basins of the 1st order, according to the characterization of hierarchical position of fluvial systems and topological ordering of hydrographic basins developed by [68]. Sub-basin 2 occurs in a micro-basin of 2nd order and sub-basin 3 on a streams junction of 3rd order (Figures 3–5). The discharge flow monitoring took place from January 2012 to December 2016. The data were collected monthly on the same point and hour sequence, from 7:00 a.m. to 11:00 a.m.

The methods used for flow measurements in watercourses are necessary to evaluate the passage of water, to obtain the measurement of the average speed of the flow in the section area and the net discharge [69–71]. The measurement of water flow in a watershed can be performed using different methods and instruments. Conventional methods include immersed current meters, which obtain the average flow velocity of the watercourse section. The velocity is then multiplied by the corresponding area, and the sum of these products will result in the average flow of the watercourse.

The measuring of flow with immersed current meters as used in natural watercourses or in artificial channels is frequent in research in Brazil. Structures are built on the riverbed, such as openings of geometric shape, as was done for this research, through which the water will flow, facilitating the measurement of the flow in the determined section [69–71]. The spillways can be: rectangular, trapezoidal, triangular, with a thin or thick threshold [69–71].

Another way of measurement could be by the artificial tracking method with the use of chemical or radioactive tracers [70,72–74]. The use of these tracers is appropriate for small watercourses, as they are considered low cost, easy to handle, presenting satisfactory results. In the past, radioactive tracers such as tritium have been used in large rivers, while tracers with fluorescents (dyes) have been commonly explored in the USA [70], especially in studies focused on groundwater [74]. However, the use of tracers introduces some problems, if the water is turbid. Suspended sediments can easily absorb some tracers and, in addition, there are normative restrictions regarding radiological protection (e.g., with the use of tritium) [70,72,73].

In this study, we opted for the electronic flow meter method (model ISCO 2150), represented as a "smart" probe". This sensor uses the concept of digital electronics. The analog level is digitized inside the sensor itself to avoid electromagnetic interference (https://www.clean.com.br/Produto/Detalhe/56). The instrument allows the monitoring of the flow at various points along the cross section of a river, in a short period, since the communication of the device with microcomputers is direct, with the transfer of the calculated flow data automatically done, substantially reducing the time necessary to fill in spreadsheets in the field and the digitization of these data. The major disadvantage of these instruments is the acquisition cost.

The discharge flow was measured with an ISCO 2150 portable flow meter (Figure 6), which allows the measurement of total flow (m3), flow speed (m/s), and level (m) of water. This equipment has doppler continuous wavelength technology to measure the average flow velocity (https://www.clean. com.br/Produto/Detalhe/56). The equipment has a cable with a sensor (flow meter) and a communicator cable for data transmission to a notebook. The flow was calculated using Isco Flowlink 5.1 software. The input parameters were configured considering the shape of each stream channel: round, U-shaped, rectangular, trapezoidal or elliptical; and by width section measurement (cm) of the watercourse (Figure 6a). The sensor was placed in each stream (Figure 6b,c) for six minutes in each sub-basin. The software automatically stored data every 30 s, and after finishing flow reading, the data were exported and saved (Figure 6d).

**Figure 6.** Flow data collection with the portable ISCO 2150 m: (**a**) width section measurement (cm) of a watercourse; (**b**) input parameters; (**c**) sensor in a watercourse; (**d**) reading flow and data storage.

In order to define the sampling points for monitoring and collecting the flow data, the area was characterized as natural sections. The 1.1 (a), 2 (c), and 3 (d) points (Figure 4) were carefully selected to determine water drainage in stream tributaries, that is, specific points where there was no separation or streambed dispersions of water flow for different parts or directions.

The flow results were subject to descriptive statistics comprising temporal and spatial analyses. The averages were showed in boxplots, with graphical representations of discharge flow distribution with a numerical summary of water flow (m3/s). For statistical treatment, the data were processed in software R, considered a tool for statistical computing and graph production (https://www.rproject.org). The pluviometric data were collected from the Integrated Center for Agrometeorological Information—CIIAGRO (http://www.ciiagro.sp.gov.br), which provides daily data by location in online system.

#### *3.2. Thematic Maps*

The land use and land cover (LULC) thematic maps of the sub-basins were elaborated using geographic information systems (GIS) and remote sensing techniques by using a visual interpretation of LULC on high-spatial-resolution orbital images.

The ESRI's ArcMap and Google Earth Pro GIS software packages were used to view, edit, create, and analyze geospatial data in hydrological and environmental studies in watersheds [26,27,72,73]. The maps' base information was compiled from spatial databases, using maps published by the

Brazilian Institute of Geography and Statistics (https://ww2.ibge.gov.br) at a scale of 1:100,000, and a digital elevation model (DEM) that was obtained from ASTER GDEN V2 satellite image.

The geoprocessing resources served as a basis for the mapping of the compositions of land use and land cover and identification of the points of collection of the water flow and of the watershed under study. The interpretative analysis of each land use was carried out based on the exposure of colors for different spectral waves of in image, with the representation and definition of the classes in regard to sub-basins' soil cover. The classification procedure was carried out based on research works developed by [27,54] and followed the characterization of environments that reflect the watershed ecosystem upstream of the water catchment point.

The location of the sub-basins, the segmentation of drainage networks and topographic dividers were carried out, and the polygons were vectorized on each LULC. The points, lines, and polygons were plotted on Google Earth Pro software, and, subsequently, the units mapped in KML format were transferred to ArcGis®, using the KML to Layer command conversion tool. The calculation of the areas for each LULC was performed in ArcMap software using the Area tool, in order to tabulate the differences of areas in percentage, in each sub-basin.

#### **4. Results and Discussion**

The headwater volumetric flow rate of water was characterized in terms of volume of fluid per unit of time (L/s) and space (sub-basins), which passed through the stream cross-section of each drainage sampling point from sub-basins: 1.1, 1.2, 2, and 3 of the Olaria stream watershed (Figure 3). Values on descriptive statistics of discharge can be seen in Table 1.

The discharge showed a range of values from 0.09 L/s (sub-basin 1.2) to 107.99 L/s (sub-basin 3), with mean values of 2.95 L/s (sub-basin 1.1); 0.92 L/s (sub-basin 1.2); 4.84 L/s (sub-basin 2); and 13.91 L/s (sub-basin 3). The standard deviation (SD) indicates the high values spread out over a wider range (1.32 to 20.23 L/s), in a high coefficient of variation (CV), positive asymmetry (γ1 > 0) and kurtosis coefficient.

The analysis of morphometric characteristics of each sub-basin provided a quantitative description of geometric aspects and slope of each sub-basin (Table 2).


**Table 1.** Descriptive statistics of discharge (L/s) of sub-basins, from 2012 to 2016.

Max—maximum value; Min—minimum value; Mean—mean value; Med—median value; SD—standard deviation; CV—coefficient of variation; AC—asymmetry coefficient; KC—kurtosis coefficient.

**Table 2.** Morphometric characteristics of sub-basins at the Olaria stream watershed.


Symbols: L—maximum length; W—maximum width; Sl—stream length; Ha—high altitude; La—low altitude; S—slope.

The data obtained from dimensional morphometry (Table 2) showed that 1st-order sub-basins (1.1 and 1.2) have similarity values in terms of area (80.9 ha and 77.8 ha), perimeter (3.54 km and 3.97 km), and stream length (0.82 and 0.80 km). The highest altitude (H) was identified in sub-basin 2 (615 m), and the lowest was on point 3 (520 m). The greatest slope (S) was on sub-basin 2 (7.4%), followed by sub-basin 1.1 (6.4%), sub-basin 1.2 (5.7%), and sub-basin 3 (5.4%). The 2nd-order sub-basin had an area of 111.2 ha and a perimeter of 4.42 km.

The headwaters (sub-basins 1.1; 1.2, and 2) form an essential link in the hydrological cycle of the Olaria stream watershed. There are several approaches to topological ordering of streams based on distance from water source (headwaters) upstream to downstream to sub-basin 3. From the highest points to the lowest points, there is a hierarchical position of 1.1 and 1.2 (first order), 2 (second order) in the river system (Table 2; Figure 3), after each stream confluence. Considering the fluvial system as a continuation of the main channel (sub-basin 3), as observed in Figure 3, the headwaters are the smaller streams, of 1st-order of magnitude [68], which are the positive integer used in geomorphology and hydrology to indicate level of branching in a river system.

The headwater upstream points of first order (sub-basins 1.1, 1.2, and 2) are the extreme tributaries on the Olaria Watershed. The Strahler order is designed to morphology of a watershed and forms the basis of important hydrographic indicators of its structure. Its base is initial discharge flow line, which characterizes spring water.

The sub-basins of 1st order of magnitude (1.1 and 1.2), with similar areas (80.9; 77.8 ha), presented different hydrological behavior regarding amplitude of volumes over time (along month and years) and space (Figure 7).

**Figure 7.** Sub-basins 1.1 and 1.2 discharge flow in the Olaria stream watershed, State of São Paulo, Brazil.

The minimum, first quartile (Q1), median, third quartile (Q3), and maximum values from sub-basins 1.1 and 1.2 showed that the water source of sub-basin 1.1 has higher volumetric flow rate of water transported through the given stream cross-sectional area than sub-basin 1.2, which showed the mostly outlier values along the monthly years, mainly in March (Month 3) (Figure 7).

The distribution, its central value, and its variability showed a different behavior of concentration of water flow on sub-basin 2, which has two stream branches and is along the drainage network (sub-basin 3) over the monthly monitored period (Figure 8).

**Figure 8.** Sub-basins 2 and 3 discharge flow in the Olaria stream watershed, State of São Paulo, Brazil.

The distribution of annual flows in sub-basins 1.1, 1.2, 2, and 3 (Figures 7 and 8) in the Olaria stream watershed indicates the values measured on an interval scale (1.1; 1.2; interval 0–8 L/s); (2; interval 0–15 L/s); (3; interval 0–100 L/s), from month 1, January, to month 12, December. Those headwaters are considered perennial, that is, they produce water throughout the year, but with flow rates varying throughout it. If two discharges from the same order stream are merged (Figure 8, sub-basin 2), the resulting discharge flow increases from 8 L/s to 15 L/s. The 2nd order of magnitude in which two headwaters from two first-order watercourses merge (sub-basin 2) showed the average value of 4.84 L/s (Table 1). Based on volume of water from confluence (the point where two 1st-order rivers merge), a minimum value of 1.44 L/s and a maximum of 17.80 L/s were observed over 5 years.

The freshwater provided from the Olaria stream on point 3 (Figure 8; sub-basin 3) showed an average flow rate of 20.47 L/s that will be used for consumption to urban areas downstream of the São Domingos River, crop irrigation systems, and multiple uses. Through the Olaria stream watershed, the discharge flow is the environmental transport on each open channel (i.e., free surface conduits of water). The water storage and control of sub-basins within headwaters (points 1.1, 1.2, and 2) determine a constant condition of water production over time.

The water discharge recognizes water source on a watershed, which allows defining the main flow and its contribution to environmental services considering the drainage networks of 1st order of magnitude. The valuation of discharge in headwaters as an ecosystem service will benefit and control the sub-basins in a watershed functioning. In this sense, discharge can be a variable that allows development of inductive and not only repressive public policies and can be considered in the importance of ecosystem services (ES) promoted by maintenance and conservation of water resources, mainly on changing the principle of "polluter pays" to "provider-receiver", in payments for environmental services (PES) schemes [53,54].

As mentioned, this program presents an importance of balancing the dynamics of habitats and ecosystems and favors maintaining, recovering, and improving environmental conditions [54,57]. So, there is a need to recognize the anthropogenic landscapes on distribution of land use and land cover (LULC) of uphill on headwaters areas. The percentage values of LULC varied in each sub-basin, with predominance for cropping and forest cover on sub-basins 1.1, 1.2, and 2 (Figure 9).

**Figure 9.** Land use land cover (LULC) at Olaria stream sub-basins, State of São Paulo, Brazil.

Agricultural land use was noticeable in all sub-basins (Figure 9). Sugarcane represented 74.8% (sub-basin 2), 56.8% (sub-basin 1.2), and 24.4% (sub-basin 1.1), annual crops were 16.9% (sub-basin 1.2) and 12.4% (sub-basin 1.1), citrus culture indicated 16.5% (sub-basin 1.1), and pasture 1.4% (sub-basin 1.2). Therefore, the sub-basins have significant anthropogenic influence. In addition to considerable crop areas, sub-basins are substantially covered with native forest, 43.6% (sub-basin 1.1), 13.8% (sub-basin 2). It is important to highlight the Agroforestry System (AS) 19.2% (sub-basin 1.2), which was implemented in 2011, aiming to restore the Permanent Preservation Area (PPA) on site and enable a better functioning of ecosystems, as well as promoting environmental services.

The characterization of water availability (upstream) allows taking into account the spatiotemporal heterogeneity of the Olaria Stream watershed, regarding physical characteristics (morphometric), LULC density, and topography. The headwaters' upstream sub-basins do not consider the crop area as a single and homogeneous unit, but regions that respond to the conservationist practices installed in an agricultural system.

The monthly average flows in sub-basins 1.1 and 1.2 (1–60 months) and the runoff coefficient (C) showed different hydrological behavior over time. The runoff coefficient (C) expresses the different relationships between the amount of runoff and the amount of precipitation received (Figures 10 and 11). The total rainfall values varied each year, (2012: 1152.1 mm; 2013: 1429.4 mm; 2014: 961.2 mm; 2015: 1367.8 mm; 2016: 1376.6 mm). The lowest value occurred in 2014 and the highest value in 2013.

**Figure 10.** Precipitation (mm) and sub-basin 1.1 and 1.2 runoff coefficient (C) from January 2012 (Month 1) to December 2016 (month 60) in the Olaria stream watershed.

**Figure 11.** Precipitation (mm) and sub-basin 2 and 3 runoff coefficient (C) from January 2012 (month 1) to December 2016 (month 60) in the Olaria stream watershed.

It is important to highlight that even in months where there was no rain or there was little rainfall (May, June, and July), the headwaters continued to produce water volume in the watershed. The sub-basins' overland surface runoff after a higher precipitation shows the water running downhill

over the landscape of the Olaria watershed (Figures 10 and 11). Stream discharges are directly influenced by seasonality, years (Table 3), and climatic variability; thus, critical periods in terms of water availability must be assessed in order to guarantee a safety margin for planning and management activities [74]. Discharge can be conceptualized as a result of hydrological processes of interaction between regional precipitations, infiltration, defluvium (Figures 10 and 11 and Table 3), and physiographic conditions of the watershed (Table 2), in which man's actions directly influence proper management of biodiversity conservation and water production [75]. When precipitation (mm) falls onto each sub-basin, it starts moving according to the laws of gravity, downhill to retain water in the streams.


**Table 3.** Sub-basins 1.1., 1.2, 2, and 3 flows (m3/s) in the Olaria stream watershed.

The water storage in the basin system has increased, and the response of the basin has become more direct as the sub-basin area increases. On sub-basin 1.2, there is an area of old gully, recovered by soil conservation practice, with the implantation of a permanent preservation area and agroforestry system, in early 2011. It is important to note that, in this place, an increase in water production can be clearly seen from September 2013 (month 47, Figure 11), a period in which soil conservation practices were consolidating along the headwaters. This result demonstrates the importance of growth of the riparian forest area, and practices of terracing, no-till, and level planting that occurred in the uphill watershed. The effect of vegetation and forests contribute to regulation of the hydrobiological cycle and greater quantity and quality of water in watersheds [52,76–79].

The water availability values of headwaters (1.1, 1.2, and 2) showed the lowest average flow in sub-basin 1.2 (0.09 L/s). Non-forested springs present several problems, such as flood events, silting of water courses by uncovered soil, in addition to low water infiltration in the hydrological system [30,80]. The lack of forest cover drastically prevents water from infiltrating the soil, thereby supplying springs, impairing proper functioning of the hydrological regime [30,81]. The headwater (1.2) is located in a large percentage of agricultural areas and an important reforestation implantation in 2011 (Figure 11). In this place, there is a clear increase in water production in 2013, confirming results of the importance of growth of the riparian forest and a better water ecosystem service promoted within a set of environmental benefits [52,55].

The results indicate the importance of a forest to increase water production and base flow (Qbf), since forest floor provides greater water infiltration in the soil (IF), intercepts (Ic) the precipitated water (P) by means of crowns, trunks and roots, which infiltrate porous soil, percolating to the deepest layers of soil, supplying the water table levels, aquifers and, consequently, springs and riverbeds, favoring the regularization of the hydrological regime. Statements like these were made by several authors, among them [69,75,82–85]. In an area with forest cover, there is less direct runoff (Qds), avoiding sediment transport, erosion, river sedimentation, flooding, and decreased loss of nutrients from soil [69,86–89].

All headwaters are perennial in an intermittent streams condition. Even in the dry season, the volume of water remained regulated, with no water outages or water crises (Table 1, Figures 7–10). In rainy months (November to April), there were no flood peaks or episodes of maximum flows, which shows the great importance of native forest in regulating discharge flow and water availability for springs, in agreement with several studies conducted [52,69,75,82–85,90,91].

The groundwater occurs within sub-basins, and the water moves slowly on sub-basin 1.2 (Table 3). The precipitation is absorbed in the watershed, where it flows and becomes part of the surface water in different ways, as shown in Figures 10 and 11. The streams flow from 1.1 and 1.2 of the watershed, join the stream on point 3, and to the São Domingos river. Thus, groundwater flows toward a stream, and the sub-basins are used as the basic hydrologic unit for both surface water and groundwater planning purposes [79–83,88,89].

Most groundwater and surface water are interconnected, but in some environments, such as sub-basin 1.2, there was not enough groundwater discharge to maintain stream flow like in 1.1. The water availability on sub-basins of first order contribute to water security assessment on the São Domingos river watershed, which contributes water flow to the Turvo/Grande river, an important Brazilian watershed.

The water security assessment framework should consider a flow discharge on a basin-scale analysis, using an indicator, mainly on first-order sub-basins, as 1.1 and 1.2. The discharge from headwaters is a dimension of water flow and can be considered as an indicator of water security. The public policies should be elaborated to consider this point of view in making the water assessment framework on the basis of headwater flow discharge. Therefore, the water availability in watersheds is a measure of how the discharge will have a bearing on water security to a better land use land cover (LULC) governance.

Throughout the period, the sub-basins functioned as an impermeable container, returning the water received by precipitation and retaining part of that water in the water storage in the dam system. The water retention capacity is influenced by several factors, among which, the forest cover, physical characteristics and geomorphological factors, topography, hydraulic works present in the basin and conservationist planting practices [26,45,47,55,67,69].

The forest cover has the function of interception, storage, and reduction of runoff [67]. Concerning sub-basins 1.1 and 1.2, morphometric characteristics at the Olaria stream watershed (Table 2) and the stream flow along the months (Figure 11), the water retention in the basin's hydrological system varied according to geomorphologic factor, as the lower slope of the terrain, soil formation, and implementation of native tree species land cover along the dams. The agroforestry system on sub-basin 1.1 was consolidated on natural vegetation cover on the basin area system along the dam, and it was observed that the flow, over time, was more stable. However, in sub-basin 1.2, in the period of growth of vegetation cover (2011 to 2013), there was a lower volume of runoff. After the consolidation of the permanent preservation area (after 2013), the volume of flow water increases gradually, with a considerable increase after 2015. According to [89], the water absorption by the root system in the growth period of the vegetation cover increased the time of permanence of water in the hydrographic basin, causing events of less volume of water flowing to the basin from the evapotranspiration processes.

With identical amounts of precipitation, the sub-basins produced varying amounts of flow, due to different physical characteristics of the hydrographic basin and vegetation cover by area. Another factor is the topography of the terrain, which may have influenced the water storage capacity on these. In areas with a higher slope (sub-basin 2; 7.6%), there was less storage capacity than flatter areas, sub-basin 1.1 (6.4%), and sub-basin 1.2 (5.7%).

The discharge flows showed space–time variation in magnitude between studied headwaters sub-basins 1.1 and 1.2, on water availability, assessed based on average net discharges (Figure 11). Another factor that can be considered is the hydraulic works present in the basin, intended to contain the runoff (dams), sub-basins 1.1 and 1.2, which result in a reduction in the maximum flow of a basin, in view of the water storage in dams. In the Olaria watershed, the lowest flow was always observed in sub-basin 1.2. It is important to note that, in this place, there are dams that were built to stabilize a gull and to recover the local spring. These dams have drainage channels to control the speed of the water in each weir.

The new Brazilian Forest Code, Law No. 12,651 of 2012 [92], currently in force, Art. 41., establishes that the Federal Executive Government will consider the "Support and incentive program for conservation of environment" for ecologically sustainable development. However, it does not declare a deadline for such action and effectiveness. In Art. 41., mentioned above, it is reported the actions that the legislation intends to apply, as well as the monetary compensations and incentives, Items I and II, namely:

I—payment or incentive for environmental services such as remuneration, monetary or not, for the activities of conservation and improvement of ecosystems and generate environmental services, such as: sequestration, conservation, maintenance and increase of stock and reduction of flow carbon, conservation of scenic beauty and biodiversity, conservation of water and water and soil services, climate regulation, cultural enhancement and traditional ecosystem knowledge, maintenance of Permanent Preservation Areas, Legal Reserves, and restricted use.

II—Compensation for the environmental conservation measures necessary for fulfillment of the objectives of this Law, namely: Obtaining agricultural credit (with lower interest rates, limits and longer terms than in the market); Hiring of agricultural insurance; Granting of tax credits (by deducting the Permanent Preservation Areas, Legal Reserve and restricted use of the calculation basis for the Tax on Rural Territorial Property—ITR); Financing lines for the recovery of degraded areas, projects for the preservation of native vegetation, protection of species of native flora threatened with extinction, forest management and sustainable agroforestry carried out on the property or rural areas; Tax exemption for the main inputs and equipment.

Therefore, Brazil is awaiting regulation from the Federal Executive Government to have a Program for Payments for Environmental Services in force in the Law. However, compensation for the generation of ecosystem services already occurs and is distributed throughout the national territory, through actions promoted in programs and projects, with federal, state, and municipal government support, or support from private-sector companies.

Some notable examples that can be mentioned are: Water Producer Program of the National Water Agency—ANA (https://www.ana.gov.br/programas-e-projetos/programa-produtor-de-agua), being executed to initiatives by city halls, basin committees, or sanitation companies. The rural producers interested in conserving springs on their property and other priority areas for water production should seek out these institutions for applications on registration in the program and financial contemplation.

Another recognized program is the Conexão Mata Atlântica project (Project for the recovery and protection of climate and biodiversity services in the southeastern corridor of the Brazilian Atlantic Forest), which benefits, through the PES financial mechanism, rural owners who adopt conservation actions native forest, recover degraded areas, and implement sustainable production practices. This program is supported by the federal government, through the Ministry of Science, Technology, Innovations and Communications (MCTIC), and by the governments of the states of Rio de Janeiro, São Paulo, and Minas Gerais. The resources are made available by the Global Fund for the Environment (Global Environmental Facility—GEF), in order to implement actions to encourage the

recovery of conservation of ecosystem services. For participation, interested parties must submit the required documentation through a public selection notice (https://conexaomataatlantica.mctic.gov.br/ cma/o-projeto/o-que-e).

There are also the Payment for Environmental Services Projects, provided for the Forest Remnants Program of the State of São Paulo—Brazil, which comply with State Law 13,798/2009 (State Policy on Climate Change) and State Decree 55,947/2010 [93,94]. These projects cover various types of related environmental services: the conservation of forest remnants; recovery of riparian forests and implantation of native vegetation for the protection of springs; planting seedlings of native species and/or implementing practices that favor natural regeneration for the formation of biodiversity corridors; reforestation with native species or with native species intercropped with exotic species for sustainable exploitation of wood and non-wood products; implementation of agroforestry and silvopastoral systems that include the planting of at least 50 individuals of native tree species per hectare; implantation of commercial forests in areas contiguous to the remnants of native vegetation; and management of forest remnants to control competing species

Other initiatives comprise conservation-planting practices. Conservation practices in agricultural crops include the use of terracing, or monitoring of level curves, to direct runoff (reducing slope) and avoid erosion and damage to crops, contribute to better conditions infiltration of water in the system, in low- or medium-intensity rains [89]. Sub-basins 1.2 have a percentage of agricultural areas that receive conservationist planting treatments, which contributes to greater water infiltration into the hydrological system. The water moves within the hydrographic basin in an infinite number of superficial and/or underground trajectories, in which there is an accumulation of groundwater and less lateral flow.

The watershed ecosystem is a fundamental unit to provide an environmental service considering natural resources. The PES programs should be applied in headwaters on sub-basins in a watershed. The scheme will bring several benefits for society, balancing the dynamics of habitats and ecosystems, and favoring maintenance, recovery, or improvement of environmental conditions [54,57,95].

The sale of a best cropping production considering best management practices and an agricultural production conducted within sustainability (e.g., agroforestry production systems) brings benefit to a collectivity. The "provider-receiver" principle is based on a landowner that will offer an environmental service which generates a benefit of a better quality and quantity of water to society. Therefore, the producer that produces in a headwater sub-basin considering agricultural practices will have the right to be remunerated for maintaining soil and land on its best quality, decreasing the erosion process and preserving the natural resources as soil and water in a better quality.

The control instrument of the "provider-receiver" principle should be a land use policy. The landowner must be motivated to include the best management practices uphill headwaters as an environmental service in their decision-making regarding a better land use land cover (LULC). The conservation of land should be a financially attractive option. The productive practices with a soil and water better management and conservation should be considered as income generation on furnishing environmental services. The PES is planning a legal framework to define the economic activity and it is an instrument of public policies at the least cost to society. The ES is determined by a market established between private and public institutions to establish a market for ES compensation.

#### **5. Conclusions**

The hydrological discharge of headwater sub-basins showed space–time variation in magnitude. The water flow discharge of forested protected headwater on 1st order of magnitude produced a better quantity of water. The Olaria stream watershed has water availability, assessed based on average net discharges (flow) and indicates a potential for environmental services payment programs, with respect to the "provider-receiver" principle, and water security schemes.

The hydrological discharge of the headwater sub-basins presented a flow range from 0.09 L/s (sub-basin 1.2) to 107.99 L/s (sub-basin 3) or from 8 to 9330 m3/day. Sub-basin 1.1 was the one that remained stable throughout the monitoring period, with an average flow of 2.95 L/s or 255 m3/day. Sub-basin 1.2 showed an increase in volume over time, ranging from 0.3 to 1.9 L/s or 25.9 to 164.1 m3/day.

Gully recovery management, reforestation with native trees, and agroforestry system in sub-basin 1.2, concomitant with the annual increase in flow during the recovery works, is an example of an activity indicated for receiving payment for the environmental service provided, discharging more water and quality for the owner downstream of the watershed.

This study contributes with an example of water security assessment in Brazilian watersheds, for the feasibility of paying for environmental services for the water-producing sub-basins.

**Author Contributions:** Conceptualization, M.B.L.S., T.C.T.P., and A.L.M.M.; methodology, M.B.L.S. and T.C.T.P.; software, L.F.S.F. and R.C.A.C.; validation, T.C.T.P., L.F.S.F., and F.A.L.P.; formal analysis, F.A.L.P., T.C.T.P., and M.Z.; investigation, M.B.L.S. and T.C.T.P.; resources, A.L.M.M. and T.C.T.P.; data curation, M.Z., L.F.S.F., and T.C.T.P.; writing—original draft preparation, M.B.L.S., M.Z., and T.C.T.P.; writing—review and editing, F.A.L.P. and M.C.L.; visualization, F.A.L.P. and M.C.L.; supervision, T.C.T.P. and F.A.L.P.; project administration, T.C.T.P. and M.B.L.S.; funding acquisition, A.L.M.M. and T.C.T.P. All authors have read and agreed 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); Agência do Ministério da Ciência, Tecnologia, Inovações e Comunicações (MCTIC); and Land Use Policy Brazilian Group (PolUS). For the author integrated in the CITAB Research Centre, the research was further financed by 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 research was further financed by National Funds of FCT—Portuguese Foundation for Science and Technology, under the project UIDB/QUI/00616/2020.

**Acknowledgments:** The authors would like to thank: the Coordination of Improvement of Higher Education Personnel (CAPES) for the scholarship; the National Council for Scientific and Technological Development (CNPq) for Research funds; the Agronomy (Soil Science) Post-Graduation Program from the School of Agricultural and Veterinary Sciences, São Paulo State University (UNESP); The Land Use Policy Research Group—PolUS. Mariana Bárbara Lopes Simedo (M.B.L.S.), Teresa Cristina Tarlé Pissarra (T.C.T.P.), Antonio Lucio Mello Martins (A.L.M.M.); Maria Conceição Lopes (M.C.L.), Renata Cristina Araújo Costa (R.C.A.C.), Marcelo Zanata (M.Z.), Luís Filipe Sanches Fernandes (L.F.S.F.), and Fernando António Leal Pacheco (F.A.L.P.).

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


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

© 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/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18