Simulations of the climate change and its effect on water resources in the Palma River basin, Brazil

With the growing water demand for various human activities, the water resources sustainability in a basin can be affected. The phenomena representation that occurs in a water basin using models is fundamental for the planning of water resources. With these models, it is possible to understand the natural processes and analyze the system response for different scenarios and allow the management to make decisions (Goes 2012). Using the system dynamics methodology, it is possible to represent complex systems and to evaluate their behavior over time, thus environments in which one needs to manage water resources can be modeled and simulated (Fan et al. 2018).

The scarcity of flow monitoring in some sites of interest, and the need to simulate system changes have led to the development of hydrological models (Caldeira et al. 2018) that are increasingly being used as a tool for water resource management and decision-making (Beskow et al. 2016), and they are suitable for evaluating different scenarios such as climate change.

Coupled atmosphere–ocean global circulation models (GCMs) are used to project possible future climate conditions, considering higher and lower greenhouse gas emission scenarios. Projections from GCMs have coarse spatial resolution for impact assessments and for regional scales, so these are refined by Regional Climate Models – RCMs (Maraun et al. 2017). Hydrological models associated with climate models can result in a future projection of flows, assisting in the management of water resources.

Based on the hypotheses that future climate change will affect water availability and demand, this study aimed to evaluate the sustainability of water resources in the Palma River basin, Tocantins, Brazil, using a system dynamics approach, throughout STELLA software combined with water supply simulated using the SWAT model, considering different scenarios of climate changes in association with different water demand in the basin, caused by the expansion of irrigated agriculture.

The study area corresponds to the Palma River basin (BHP), presented in Figure 1, which is in the Southeast region of the Tocantins state (TO) and belongs to the Tocantins River Hydrographic System (right bank). The Palma River basin is located between the South parallels 13°03′ and 11°24′ and the West longitude meridians 47° 53′ and 46° 04′.

The Palma River basin has a drainage area of 17,468 km², and for its delimitation, it was taken as a control section the Barra do Palma fluviometric station, which is located on the Palma River, Paranã city – TO (Rodrigues 2017).

To calculate the basin water demands, data were obtained from the population and agricultural censuses of the Brazilian Institute of Geography and Statistics (IBGE), from the Population Projections (Brazilian Institute of Geography and Statistics 2020), Municipal Livestock Research (PPM) (Brazilian Institute of Geography and Statistics 2016a) and Municipal Agricultural Production (PAM) (Brazilian Institute of Geography and Statistics 2016b) platforms. The data concerning to the growth of irrigation were obtained from information made available by the Brazilian Association of Machinery and Equipment Industry (ABIMAQ) (Brazilian Machinery and Equipment Industry Association 2018).

In a GIS software, having the delimitation of the Palma River basin and the municipalities located in the basin, a clipping was made obtaining the area of each municipality within the basin. Then having the number of inhabitants, animals, and irrigated land in each municipality, this number was estimated for the area of each municipality within the basin, and the sum of all municipalities was made, obtaining finally the total number of inhabitants, animals, and irrigated land area in the basin.

The average water consumption per inhabitant was obtained from the National Water and Sanitation Agency – ANA (2012), and water consumption by animals was obtained from Matos (2003). The water requirements for crop irrigation were simulated using the AquaCrop software (Steduto 2009, 2012).

The flowchart in Figure 2 presented the steps that were followed in the simulation of the water resources in the Palma River basin. Rodrigues (2017) applied the climate projections from the Eta/HadGEM2-ES and Eta/MIROC5 models, for the RCPs 4.5 and 8.5 emission scenarios, in the SWAT model to simulate runoff and potential evapotranspiration in the basin. The potential evapotranspiration data were used as input in the AquaCrop model, enabled the generation of future water demand data for crop irrigation. Finally, the irrigation, livestock, population, and environmental demands, along with the runoff simulated by SWAT were used as input data in the STELLA software for the simulation of water supply and demand in the basin under each climatic change scenario.

The future runoff and potential evapotranspiration data used in this study were simulated using the SWAT model by Rodrigues (2017), who evaluated in his study the anthropogenic impacts on the hydrological regime of tributaries of the Tocantins River, including the Palma River basin.

In Rodrigues' (2017) work, the impacts of climate change on runoff in the BHP were estimated from hydrological simulations based on climate projections from the Eta-CPTEC/MIROC5 and Eta-CPTEC/HadGEM2-ES models. These projections were produced for the present climate (1961–2005) and for RCPs 4.5 and 8.5, over the 21st century (2007–2040, 2041–2070 and 2071–2099) at 20 km resolution.

To simulate different crops water requirements for irrigation was used version 6.1 of the AquaCrop software (Food and Agriculture Organization of the United Nations 2018), developed by the Food and Agriculture Organization of the United Nations – FAO (Steduto 2009, 2012).

The input data referring to climate in AquaCrop were the potential evapotranspiration obtained in simulations performed by Rodrigues (2017) through the SWAT model, using the GCMs Eta-CPTEC/MIROC5 and Eta-CPTEC/HadGEM2-ES, and the RCPs 4.5 and 8.5 emission scenarios for each model.

A standard soil type presented by the software was chosen, consisting of a medium-textured soil, with great depth and no water input from the water table.

The most expressive crops in the basin were analyzed according to the last IBGE census and were used in the simulations those that most use irrigation for their management such as, sugarcane, beans, corn, sorghum, and soybeans. Historical series from the platform Municipal Agricultural Production (PAM) of IBGE was used to evaluate the proportion of each crop considered in the municipalities of the basin over the years. Given the difficulty of allocating crops within the study area, it was decided to disregard the effect that this allocation would have on the results, thus considering that the crops are randomly distributed within the basin. It was chosen to separate the crops into annual and perennial crops, with sugarcane being the only representative of perennial crops. The use of sugarcane separately was due to the importance it has been gaining in the region, and for its unique behavior regarding to water requirement for irrigation in relation to the other crops studied.

The simulations were run and for each simulated crop cycle the net irrigation requirement was calculated, which represents the water requirement to irrigate the crop without considering the efficiency of the irrigation method used. For the simulation in STELLA a fixed value of 0.85 efficiency was used, this being an average value considering different irrigation systems. The gross demand was obtained from the ratio between the net demand and the irrigation efficiency.

This study uses the STELLA software version 9.0 (STELLA 2001), which uses the system dynamics methodology, and works with physical, economic, and social variables. The construction of the model has been structured through the stages of conceptualization; formalization; and simulation, which includes evaluation and exploration.

In the conceptualization stage, the planning horizon for the simulations was chosen from 2007 to 2099, with the time unit being the month. The year 2007 was chosen for the beginning of the simulations because from this year on there is no flow data available at the fluviometric station located in the basin.

Water demands in the basin comprised: domestic (urban and rural), agricultural, livestock, and environmental demands. Environmental demand is that which is necessary for the protection of the environment and ecosystems. In this study, the environmental demand was assumed to be the minimum average monthly flow with 90% probability of occurrence (Q90), observed in a series of flows from 1961 to 2005.

To calculate the population demand, the population growth rates projected by the IBGE for the state of Tocantins throughout the century were considered. And to calculate the agricultural demand, the growth rates of irrigated agriculture and livestock observed in the data made available by the Brazilian Association of Machinery and Equipment Industry (ABIMAQ) and the Municipal Livestock Research (PPM) respectively, were considered.

After the conceptualization stage, there is the formalization stage of the model using a mathematical language. For this, a stock and flow diagram (Figure 3) was elaborated, which describes the system's operation from which the model's equations are written.

If the SI is higher than 0.2, there is no or little water stress in the basin, meaning that demand consumes less than 80% of supply; if it is less than 0.2 there is vulnerability, meaning that demand consumes more than 80% of supply; and if the SI is equal to zero the supply is not sufficient to respond demand (Xu et al. 2002).

After creating the stock and flow diagram, the data were entered, such as those already simulated for the supply and those for the demands as well as the growth rates. Then the simulation duration and the time unit were established. Finally, the simulations were run, obtaining the various components behavior over time.

The STELLA software allows through its analysis tools checking of syntax errors in the equations and the analysis of consistency between the units used in the model. As mathematical relationships are established they are automatically checked and if any error is detected the equation cannot be established (Ortega & Ortega 2017). The Check Units tool was used in each variable established and also in the model as a whole in order to check the consistency between the units.

To evaluate the dimensional consistency of the model, the equivalence of the equation Supply=Demand+Stock was checked, as was done by Orellana González et al. (2008).

The validation was also done by comparing the behavior of the model's variables, such as population and agricultural demands. For example, as the population grows so does the demand for food and consequently the irrigated land. Thus, the population and agricultural demands for water are expected to increase, which was verified.

To validate the model behavior, the population and agricultural growths were compared with existing data for the period 2006 to 2016. The population data were compared with IBGE projections for the same period. And the data referring to livestock were compared with data from the Municipal Livestock Research (PPM) of the IBGE. The comparison of simulated and observed data was made by calculating the Pearson correlation coefficient.

Four different scenarios of irrigated area expansion were evaluated. The evaluated basin has been prominent, and it is presenting great potential for its development in the coming years.

The average growth rate of the irrigated area in Brazil in recent years (3.8% per year) was considered, it was obtained from data made available by ABIMAQ. It was also assigned an irrigated area growth rate of 8.5% per year, where the irrigated area reaches 50% of the total basin area, which according to the Secretary of Environment and Water Resources – SEMARH is the percentage of basin land suitable for irrigation. An intermediate irrigated area growth rate of 6% per year and a scenario with no growth in the irrigated areas was also considered.

For each irrigated area growth rate, the SI was simulated considering the runoff generated in SWAT by Rodrigues (2017) based on the climate projections of the Eta/HadGEM2-ES and Eta/MIROC5 models in RCPs 4.5 and 8.5 emission scenarios. Thus, a total of 16 scenarios was evaluated.

Correlation tests were performed, using the XLSTAT software in its test version, to detect correlations between the variables of the SI in the different scenarios evaluated.

The correlation tests were used to calculate different types of correlation coefficients, between two or more variables, and verify whether the correlation is significant or not. In this study, the Pearson correlation coefficient was chosen, which is suitable for continuous data and ranges from 1 to −1. p-values were calculated, which allows testing the null hypothesis that the coefficients are significantly different from zero (XLSTAT 2020).

The XLSTAT also performs an AHC, which is used to form homogeneous groups of objects based on their description by a set of variables, or from a matrix that describes the similarity between the objects (XLSTAT 2020). One of the results generated was the dendrogram, which showed the progressive grouping of the data, so one could get an idea of the number of classes in which the data can be grouped.

To calculate the dissimilarity between the groups, Ward's method (1963) was used, this method aggregates two groups so that the inertia within the group increases as little as possible to keep the groups homogeneous. This method can only be used in cases with quadratic distances such as the Euclidean distance, which was used in this study.

To check for trends in the simulated SI series, the Mann–Kendall test (Mann 1945; Kendall 1975), which is a nonparametric test for trends in time series, was used. The null hypothesis H0 for this test is that there is no trend in the series. The tests are based on the calculation of Kendall's tau correlation coefficient, which can be interpreted as the difference between the probabilities that the variables vary in the same direction and the probabilities that the variables vary in the opposite direction (XLSTAT 2020).

A significance level of 5% was considered, so the hypothesis of no trend was rejected when the p-value was less than the significance level.

The Mann–Kendall test requires independent observations, meaning that the autocorrelation should not be significant, as the variance of the S-statistic may be underestimated (XLSTAT 2020). To check for autocorrelation in the series, the test of Yue & Wang (2004) was chosen, which performs best when there is a trend and autocorrelation. In this method the Sen's slope estimator (Sen 1968) is calculated, the closer this is to 0, the lower the trend is, and the sign indicates whether the trend is increasing or decreasing.

Comparing data simulated from STELLA, with the observed data from IBGE, a Pearson correlation coefficient of 0.8 was found for livestock, and 0.95 for the population. According to Vieira (2011), if the correlation coefficient is between 0.75 and 1.00 the correlation between the variables is strong, which was verified.

The behavior of precipitation over time, for the different climate changes scenarios, simulated by Rodrigues (2017) is presented in Figure 4.

Table 1 shows the results of the Mann–Kendall trend test for the precipitation series for each climate change scenario.

Only in the HadGEM 4.5 scenario was no trend observed in the precipitation series, in the other scenarios there was a tendency for precipitation reduction throughout the time. It can be observed that the trend is even greater in MIROC 4.5 and 8.5 scenarios, as the p-value is very close to zero. Figure 4 also show that the angular coefficient of the trend line in MIROC 4.5 and 8.5 scenarios are larger, indicating a greater reduction in precipitation over time.

To check for the existence of autocorrelation, the calculated Sen's slope estimator showed a negative value, which indicates that the trend is downward, and is higher in the MIROC 4.5 and 8.5 scenarios.

The behavior of the irrigated area obtained from the STELLA software for the different growth rates over the simulation time is shown in Figure 5.

When considering the growth rate of 3.8% per year, the total irrigated area reaches approximately 12,000 hectares at the end of the simulation. Considering the intermediate growth rate of 6%, the irrigated area reaches almost 90,000 hectares in 2099, and considering the growth rate of 8.5%, the irrigated area reaches approximately 900,000 hectares at the end of the simulation.

Figure 6 shows the SI in each climate change scenario for the irrigated area growth rates of 3.8%, 6%, and 8.5% per year, and without irrigated area expansion. The line of SI=0.2 was represented because if this is less than 0.2 there is a situation of vulnerability of water resources in the basin, where the demand consumes more than 80% of the supply.

Table 2 presents the summary of the results for the SI in the different climate change scenarios, for each expansion of the irrigated area simulated.

As observed in all scenarios of irrigated area expansion, HadGEM 8.5 presented a lower SI mean, a higher number of events with SI less than 0.2, a longer consecutive period with SI <0.2, and a shorter consecutive period with SI >0.2, which indicates this scenario presents higher impacts on the basin water resources than the others.

Comparing the irrigated area expansion scenarios, it a reduction was observed in the mean and median SI as the growth rate of the irrigated area increases, as well as an increase in the number of IS <0.2 events, indicating that the expansion of the irrigated area can affect the basin's water resources, leading to its unsustainability.

Table 3 presents the Pearson correlation matrix for the SI in the different climate change and irrigated area expansion scenarios.

In the correlation matrix, the most significant correlations are indicated in bold, with those closest to 1 being the strongest positive correlations, and those closest to −1, the strongest negative correlations. Note that there was a stronger positive correlation between the variables: SD and SI mean; SI cv and SI events <0.2; and SI >0.2 longest period and SD. A strong positive correlation was also observed between the SI >0.2 longest period and SI mean, indicating that the higher the SI mean, the longer the consecutive period when the condition in the basin is sustainable.

There was a greater negative correlation between the variables: SI cv and SI mean; and SI <0.2 and SI mean, indicating that the lower the SI mean, the greater the number of unsustainable events, and vice versa. A strong negative correlation was also observed between the number of events with SI <0.2 and the longest period of SI >0.2, indicating that the greater the number of unsustainable months, the shorter the consecutive period with sustainable SI.

Figure 7 shows the dendrogram with the scenario groups formed from the correlation matrix.

It can be observed that two larger groups have been formed, with one comprising the HadGEM 8.5 scenarios, and the other group for the other ones. The HadGEM 8.5 scenarios are the most pessimistic, presenting the lowest SI mean value, as well as the highest number of events with SI <0.2, thus differing from the other scenarios.

When analyzing the other group formed, the dissimilarities were very small, so much so that in the extension of the dendrogram the scenarios in this group are all equal. The HadGEM 4.5 scenarios presented a small dissimilarity from the other scenarios, as the HadGEM 4.5 scenarios showed higher SI mean value, lower number of unsustainable events, and longer consecutive period with sustainable events. The MIROC 4.5 and MIROC 8.5 scenarios were very close, as the results found in these scenarios were similar.

Table 4 shows the results obtained in the Mann–Kendall trend test for the simulated SI series, considering the different climate change and irrigated area expansion scenarios.

According to the trend tests, in all irrigated area expansion scenarios evaluated, only the HadGEM 4.5 climate scenario showed no trend in the simulated SI series, since the p-value was greater than the significance level of 0.05, so the hypothesis H0 could not be rejected. In the other scenarios, the p-value was less than the significance level, thus the hypothesis H0 was rejected, which indicates that there is a trend in the series. Thus, there is a tendency for SI to decrease over time in most of the climate scenarios evaluated.

Sen's slope estimator, used to verify the existence of autocorrelation in the series, was very close to zero in all scenarios, which shows that the trend is very small.

The percentage of unsustainable months in which the SI was less than 0.2 was calculated for all years in all scenarios. Figure 8 shows the percentages of unsustainable months for each scenario over the simulation time, considering the different growth rates of irrigated area.

It is observed that, although the HadGEM 8.5 scenario was more unsustainable, regarding the SI mean and the number of events with SI <0.2, the MIROC 4.5 and 8.5 scenarios presented higher angular coefficient of the line, regarding the percentage of unsustainable months. Meaning, in these scenarios, there was greater growth in the number of months with SI <0.2 until the end of the simulation, which can be attributed to the fact that there was a greater reduction in precipitation in these scenarios, which was observed by Rodrigues (2017).

Table 5 present the results of the Mann–Kendall trend test, for the series of percentages of unsustainable months over the simulation years, for the different growth rates of irrigated area.

According to the trend tests, in all the irrigated area expansion scenarios evaluated, only the HadGEM 4.5 climate scenario showed no trend in the simulated series. In the other scenarios, there was an increasing trend in the percentage of unsustainable months over the years.

About the autocorrelation in the series, the calculated Sen's slope estimator was positive, which indicates that the trend is increasing and is higher in the MIROC 8.5 scenario.

In the Rodrigues (2017) work, which simulated the runoff data used in the present study, the projections of the Eta/HadGEM2-ES and Eta/MIROC5 models for RCPs 4.5 and 8.5 in the BHP showed trends of runoff reduction, and the most severe reduction may occur in the period 2071–2099. In the present study, a trend of reduced SI and an increase in the percentage of unsustainable months was observed in most of the scenarios evaluated throughout the simulation, which may have been caused by the reduction in the runoff, as well as by the increase in water demand.

Also, according to the same author, the Eta/HadGEM2-ES model for RCP 8.5, was the one that presented the highest projected hydrological impacts for the studied basin, which was verified in the present study, as Eta/HadGEM2-ES for RCP 8.5 was the model that presented the lowest SI mean value and the highest number of events with SI <0.2, in all simulated irrigated area growth rates.

Souza et al. (2010) analyzed the hydrologic behavior in the Entre Ribeiros stream basin using the STELLA software. In the climate change scenario, which comprises future increases in precipitation and temperature, in the year 2039, the SI was 0.2, indicating that there would already be hydric stress in the system. In the present study, an SI mean of 0.43, 0.08, 0.27, and 0.16 for the HadGEM 4.5 and 8.5, and MIROC 4.5 and 8.5 scenarios, respectively, was found for the BHP considering the current irrigated area growth rate of 3.8% per year in the year 2039. This indicates that the water resource situation in the basin would also already be vulnerable in this year for the 8.5 emission scenarios, which are the most pessimistic.

Xu et al. (2002) simulated the SI using STELLA, for a sub-region of the Yellow River basin in China. In a scenario of a 20% increase in irrigated area, the SI was 0.29, 0.18, and 0.11 for the years 2010, 2020, and 2030, respectively, denoting that the sub-region is already experiencing water scarcity. For the BHP considering the irrigated area growth rate of 6% per year, it was found in the year 2010 the SI mean of 0.52, 0.01, 0.35, and 0.45; in 2020, the SI mean of 0.11, 0, 0.43, and 0.54; and in 2030, the SI mean of 0.08, 0.18, 0.16 and 0.14, for the HadGEM 4.5 and 8.5, and MIROC 4.5 and 8.5 scenarios, respectively. These results present that in the coming years the BHP may also present water scarcity if there is an increase in the irrigated area, under the climate scenarios considered.

Orellana González et al. (2008), in a study in the municipality of São Miguel do Anta (MG), evaluated the water resources sustainability based on the SI. Considering an integrated scenario of climate change and an increase in irrigated area, where 50% of agricultural land is irrigated, in the years 2011, 2019, 2027, and 2035, the SI mean was 0.34, 0.32, 0.30, and 0.27, respectively. For the BHP, considering the irrigated area increment rate of 8.5% per year, where 50% of the basin lands are irrigated, the SI mean of 0.53, 0.11, 0.48, and 0.48 were found in the year 2011; in 2019, the SI mean of 0.11, 0.13, 0.41 and 0.51; in 2027, the SI mean of 0.16, 0.02, 0.37 and 0.3; and in 2035, the SI mean of 0.18, 0.09, 0.27 and 0.08, for the HadGEM 4. 5 and 8.5, and MIROC 4.5 and 8.5, respectively. Both studies showed a trend of decreasing SI over the years under climate change scenarios and increasing irrigated area, which may lead to the compromise of water resources in the studied basins in the coming years.

With the simulations of water supply and demand in the Palma River basin, considering the different climatic change and irrigated area expansion scenarios, it was possible to evaluate, based on the SI, the situation of the water resources in the basin until the end of the century.

With regard to the water supply, this presents a reduction throughout the simulation. Conversely, water demand increases, especially in the scenarios where there is a bigger expansion in the irrigated area, and in those that are more pessimistic in terms of greenhouse gas emissions. Therefore, events occur in which demand exceeds supply, especially at the end of the simulation, which causes the SI to be lower than 0.2, indicating the unsustainability of water resources in the basin.

The results found show future problems may occur in water availability in the basin, and thus can assist in the management of water resources respect to water allocation, considering the water supply to current and future users and meeting environmental demands, especially in periods of reduced runoff. The results can contribute to decision-makers, so the use of water resources in the basin can continue without affecting its ecosystem and economic sustainability.

This study was supported by the Fundação de Amparo à Pesquisa de Minas Gerais (FAPEMIG), through the scholarship for the first author fund number 5.316/15.

All relevant data are included in the paper or its Supplementary Information.