The effect of climate change on surface and groundwater resources using WEAP-MODFLOW models

In addressing management scenarios and climate changes, it is necessary to consider surface water and groundwater resources as an integrated system. In this context, the present research first simulates and evaluates surface water and groundwater simultaneously; then, it examines the possible effects of climate change on these water resources in the study area (Mahabad, Northwest of Iran). In the first stage, the WEAP-MODFLOW model was applied to a 10-year period (2006–2015) in order to take into account the interactions between surface water and groundwater and calibrate the amount of recharge and drainage from the aquifer. In the second stage, in order to study the effect of climate change on surface water and groundwater resources, we compared the micro-scale model outputs under the RCP4.5 scenario for different climate change models in the period 2021–2045. The results show that root-mean-square error (RMSE) and mean absolute error (MAE) scores are equal to 0.89 and 0.79 in unsteady conditions, respectively, which confirm the efficient performance of groundwater simulation. In addition, the results of the WEAP model based on MARE assessment criteria for calibration and validation modes are equal to 0.54 and 54.0, respectively. This finding provides evidence for the efficient performance of the simulation model. Once the interactions between groundwater and surface water were specified, the results R2 and NS suggested that indices were equal to 0.62 and 0.59, respectively, for Mahabad hydrometric station. The efficient performance of the proposed model for runoff simulation was therefore confirmed. Owing to climate change in the study period, groundwater decreased by about 1.6–1.9 m. Moreover, the amount of runoff declined from 0.1 to 0.001 MCM/month in all months except December. Unless appropriate decisions are taken to improve groundwater and strategies are applied to reduce the effect of climate change, under the present conditions this region will suffer irreparable damages in the future.


Introduction
Proper management of surface water resources (including dams and rivers) and groundwater resources requires paying special attention to all parameters of hydrologic balance. Water consumption management is vital in arid and semiarid regions. The agricultural sector is the largest consumer of water, accounting for about 70% of the consumption of fresh water resources (FAO 2017).
Water shortage and restriction of water resources is currently a major problem in Iran. Therefore, the role of water resource and consumption management is especially important due to an increasing demand for water. In addition, Iran is located in an arid and semiarid region. Population growth, urbanization and expansion of urban areas, and industrial and agricultural developments have led to an increased water demand (Hashemi et al. 2018;Ostad-Ali-Askari 2022). Under these circumstances, water resource management is crucial in order to maintain the balance between water demand and consumption . Therefore, it is necessary to have a sound understanding of the natural behavior of the hydrological system in order to manage hydrological phenomena. Since surface water and groundwater resources are the two main systems for meeting water demand, it is necessary to take their capacities and limitations into consideration (Milan et al. 2023).
Groundwater is one of the resources of irrigation water (Barati 2020a, b). Reasonable and balanced exploitation of these resources ensures their sustainable use, and their excessive exploitation ends in their depletion (Milan et al. 2023;Kayhomayoon et al. 2022;Nematollahi and sanayei 2022;Maghrebi et al. 2023;Ostad-Ali-Askari et al. 2022). Indeed, excessive extraction of groundwater affects the quality and quantity and other areas such as the environment, agriculture, and industry. It is therefore necessary to maintain a balance in using these resources in order to preserve and restore water ecosystems such as rivers and wetlands.
In recent years, ignoring the potential capacities of surface water and groundwater resources in different regions of Iran has caused various problems such as water shortage in the urban, agricultural, and industrial sectors as well as environmental damages such as loss of wetlands (Milan et al. 2023;Tork et al. 2022;Du et al. 2022). Integrated water resource management with an emphasis on conjunctive water use from surface water and groundwater resources is part of the agenda of the organizations responsible for water exploitation. Conjunctive water use implies the exploitation of both surface water and groundwater resources in order to supply water demand. Development of groundwater exploitation has considerable advantages and fewer disadvantages compared to surface water exploitation strategies such as construction of dams (Scanlon et al. 2023;Wu et al. 2022). It is more cost-effective, it does not cause problems such as deposition and precipitation, there are fewer issues related to quality, and it does not involve social and cultural challenges. The water evaluation and planning system (WEAP) model is able to simultaneously take into consideration different important factors that impact natural and man-made water resource management systems and river basin management. This model is easily accessible and is widely used on a global scale for solving similar problems (Shafa et al. 2023; Ostad-Ali-Askari 2022). These qualities have attracted researchers' attention to the WEAP model. In regions where aquifers and rivers are used in an integrated fashion, groundwater is considered as a good option for exploitation in case of a river drought. Since there is a mutual effect between the hydrological cycles of surface water and groundwater resources, water-level decline in one of the resources inevitably impacts the other (Dehghanipour et al. 2019). Consequently, surface and groundwater resources must be used and considered simultaneously in water consumption management in order to optimally exploit all the potentials of the existing water resources while preserving them at the same time. Therefore, using numerical models such as MODFLOW is essential in understanding the behavior of groundwater resources. The objective of the mathematical model of an aquifer is to simulate the natural behavior of the aquifer through a series of mathematical relations. In this simulation, the river is considered as a surface water resource and the model explores the interactions between the two water resources. If we simulate an integrated system of surface water and groundwater resources and tune it with natural conditions, we can easily study the impact of exploitation on the aquifer by changing the location, time, and amount of exploitation. Conjunctive use of surface water and groundwater resources can help preserve the existing water resources, improve the effectiveness of water resource management, and minimize the negative effects of the conventional approach, which is focused on a single water resource. Integrated use of water resources aims to increase possible exploitation amount while maintaining sustainable use of the existing water resources. This strategy reduces the costs of exploitation, makes sustainable use possible while groundwater resources are used as supplementary sources, and helps maintain surface water resources approximately within the same level at times of drought and wet years (Banihabib et al. 2017;Dehghanipour et al. 2020;Dehghanpour et al. 2019).
Climate changes affect the duration, intensity, type, and time of precipitation in different regions of the world, which in turn cause changes in water supply (Piao et al. 2010;Wheeler and von Braun 2013). Consequently, the change in groundwater supplies incites tensions in water management (Neto et al. 2021). Climate changes have a direct impact on the amount of surface discharge and an indirect impact on the storage of aquifers. There is a close relationship between the components of the water cycle and the climate system, causing them to have a direct effect on each other. The amount of runoff, river discharge, groundwater, and the intensity of flood and drought are all influenced by the most important climatic factors (i.e., temperature and precipitation). Investigating climate change and its effects on water resources in general and runoff in particular can pave the way for adopting efficient water resource management policies in the future. In the last decade, many studies related to climate change have focused on surface water and groundwater resources (Ficklin et al. 2009;Abbaspour et al. 2009;Neto et al. 2021;Azizi et al. 2021;and Zhang et al. 2012). These studies have dealt with the possible distinct effect of climate change on surface water and groundwater resources and have rarely addressed its concurrent effects on surface water and groundwater resources (Maghrebi et al. 2021). Using GMS-MODFLOW, Moghadam et al. (2022) explored the effect of climate change on surface water and groundwater resources.
Given the climatic changes during the last few decades and the effects of this phenomenon on climatic parameters, the future cyclic pattern of hydrological phenomena will not be the same as in the past. Hence, employing different climate models that can predict climate changes, integrating these models, and finally using the output of these models in hydrological models can provide a useful instrument for predicting hydrological changes in the future. The disadvantage of the conventional approach, which is focused on a single water resource, is that it fails to discover the exact impact of climate change on water resources. Indeed, surface water and groundwater are strikingly interconnected and their analysis requires an integrated perspective (Mirani Moghadam et al. 2021). In this regard, the present study aims to pinpoint the effect of climate change on surface water and groundwater resources by using the WEAP-MODFLOW model. To this end, this integrated approach is applied to Mahabad aquifer, northwestern Iran.
Mahabad aquifer has been used as a supplementary water resource for the region. In the past, Mahabad Dam has served as main surface water resources in this region. However, in recent years it supplies a smaller share of the demand due to rainfall shortage. Consequently, the exploitation of groundwater resources has intensified. A possible scenario for increasing the balance of Lake Urmia involves exploitation of surface water resources of Mahabad. In this scenario, groundwater resources must supply the water demand. Consequently, it is necessary to explore alterations in groundwater balance under unsteady conditions such as climate changes.

Study area and data
The catchment area of Mahabad River is in the south of Lake Urmia, northwestern Iran. This basin extends over 53.1524 km 2 , which comprises about 3% of the catchment area of Lake Urmia (Fig. 1). The catchment area of this river is located between 45 degrees, 25 min, and 9 s to 46 degrees, 45 degrees, and 51 s east longitude and 36 degrees, 23 min, and 51 s north latitude. Mahabad River is formed by the confluence of Bitas tributary in the east and Koter tributary in the west. Satellite images show this river is highly meandering, which is less common in other rivers flowing into Lake Urmia basin. In 1967, an embankment dam was built on this river at the junction of the two tributaries Bitas and Koter. At the same time, another dam was constructed downstream of Mahabad Dam to divert water to irrigation canals and to control flood, which happens once a year in the lower reaches of the basin. Mahabad Dam is one of the ten largest dams in Iran. Having a core built of pebble and clay, the dam has a height of 47.5 m and a length of 700 m. This dam was built on the Mahabad Chai River, and its annual water input is 339 million cubic meters. Located one kilometer from the city of Mahabad, the reservoir of this dam (extending 360 hectares) is a permanent wetland that supplies municipal and agricultural water to the city and surrounding villages (Asalpishe and Manaf Far 2017).
We have considered monthly data for surface water and groundwater modeling because it was more easily accessible and available compared with daily and annual data.

Methodology
As shown in Fig. 2, this research was carried out in four phases. In the first phase, groundwater was simulated using MODFLOW. In this phase, the quantitative status of groundwater was assessed in two modes: steady and unsteady. The purpose in this phase was to study the impact of climate change and management scenarios in the future. In the second phase, the geographical area under study was simulated using WEAP and the amounts of water shortage in each month were specified. In the third phase, using the selected model, we examined changes in precipitation and temperature for the 2025-2045 period. Downscaling was conducted using the Delta method. Eventually, the impact of climate change was applied to surface and groundwater resources. In the fourth phase, management scenarios were applied to better meet the requirements of surface and groundwater resources.

Water evaluation and planning system (WEAP)
Using the WEAP method, we modeled the water resources and uses of the downstream area of Mahabad Dam. A number of scenarios for water consumption management were also implemented in this model. In the next step, the proposed model was calibrated and validated on the basis of data from the hydrometric station of Mahabad River. The model was developed for a period of 15 years: the first 10 years were considered for calibration and the last 5 years were intended for validation.
Besides, the mean relative absolute error (MARE) was used in order to estimate the calibration and validation error in the simulated model.
where Q ci is the discharge calculated by the model in the i-th year, Q oi is the observed discharge of the hydrometer station in the i-th year, and n is the number of data.
In the upstream PI decentralized control system, the control variable (i.e., the measured discharge passing through the regulatory structure) is calculated in accordance with the error rate of water level relative to the target value; finally, the degree of the valve opening is determined. The discharge passing the structure located in the canal's headwaters enters the canal through a feedback controller and in proportion to the water demand values for each catchment area (or based on the specific water right of each area). The proportional-integral control algorithm could be applied to measure changes in the discharge passing through the regulatory structures using Eq. 2 (van Overloop, 2006;Carpenter et al. 2022).
where ∆Q (k) is the value of the controlled discharge passing through the regulatory structure in cubic meters per second in the current time step; e is the difference between the [actual] water level and the target level; k and k + 1, respectively, denote the current time step and the previous time step; and k p and k i , respectively, denote the proportional coefficient and the integral coefficient. To design the PI controller, we calculated the proportional and integral coefficients (k p and k i ) using the formula proposed by Schuurmans (Schuurmans 1997).
where A s signifies the reserve level under the exploitation conditions of the control time step, k c is the maximum resonance in the minimum flow, and R p is the maximum resonance frequency in the minimum flow of the canal range and can be measured using Eqs. 3 and 4 (Schuurmans 1997). The MPC controller predicts the control system by integrating two factors: time horizon and mathematical model. In controlling the water system through the MPC method, the state-space model is used to express the internal model. Consequently, it will be possible to compress the multivariate formulation of linear models. The objective function used for the canal system is to define: where J is the objective function and must be minimized; X refers to state variables; U refers to control actions; Q is the weight matrix for state variables; and R is the weight matrix for control actions. The optimization algorithm used in the MPC controller can limit the solution space by applying certain constraints. The constraints of this research include the amount of maneuvering of the controlled valve in each time step (u lim ) and the minimum amount of required water (x lim ). Equations 6 and 7 show hard constraints for the values of the controlled state and variable: where E and F represent the selected matrix with the initial values of 1 and − 1, so that in order to equalize the calculations, the form of the constraints will always be smaller unequal and equal.
As shown in Fig. 3, the simulation was applied to a 10-year period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) in the present research. In the next step, we included data related to Mahabad River and its tributaries, Mahabad Dam reservoir, and the amount of inflow to the dam reservoir in the model. Next, we prioritized domestic, environmental, industrial, and agricultural water consumption in that order and included them in the model separately. We also considered exploitation rate of surface water and groundwater supplies as well as the amount of return flows to groundwater and surface water resources in the model. These data were obtained based on meteorological information of Mahabad synoptic station. In this case, the amount of surface water and aquifer recharge from domestic, industrial, and agricultural sectors were equal to 0.75, 0.65, and 30 percent, respectively. Eventually, water resources and water consumption of the study area (traditional agricultural practices, industrial water demand, domestic water consumption, and environmental requirements) were simulated using WEAP model (Fig. 4). In order to estimate the calibration and validation error values in the simulation model, we used Nash-Sutcliffe coefficient of efficiency (NSE) and mean absolute relative error (MARE) (Eq. 4).
After the data were entered in the model and the river flow was simulated, WEAP resolved the water allocation problem. The objective function in WEAP model is maximization of coverage percentage of water demand of each node, which is obtained through linear planning.

Groundwater simulation
One must develop a quantitative model of groundwater in order to implement management scenarios and assess the effect of climate change along with the effect of integrated exploitation of water resources. In this regard, we drew on hydraulic, hydrological, and geological information to perform numerical simulation of Mahabad aquifer using the MODFLOW code (Harbaugh et al. 2000) and the groundwater modeling system (GMS) software. Next, we examined the exploitation status of the aquifer. One of the most wellknown quantitative tools for groundwater modeling is the MODFLOW code, which has been developed by the US Geological Survey. It is capable of modeling groundwater in steady and unsteady circumstances. To this end, we must first build a conceptual model. This model serves the following functions: determining the boundary of the aquifer and the location of exploitation wells and observation wells; determining the location of the input and output flows of the aquifer; estimating the initial hydraulic conductivity values; and determining the amount of aquifer inflows from surface and return flows. There are around 1,715 wells in Mahabad aquifer, which extract 16.23 million cubic meters of groundwater every year. The locations of these wells are shown in Fig. 3. There are also 22 observation wells (Fig. 3) to control the level of groundwater in the aquifer. Mahabad River is also included as a coverage in the conceptual model. The period from September 23 to October 22, 2010, was considered steady due to the least changes occurring in groundwater in this 30-day period. The level of groundwater for the observation wells, the discharge from the exploitation wells, and the amount of recharge from precipitation and floods were quantified under different coverages in this period. The time from October 23, 2019, to September 22, 2012, was considered as the calibration period. Also, a grid cell size of 250 × 250 was set for the aquifer.

Climate change simulation
The models of the IPCC-Fifth Climate Change Report were used to simulate climate change in the study area. Three scenarios RCP2.6, RCP4.5, and RCP8 were used to generate future (2026-2045) precipitation and temperature data. For this purpose, 14 different models of CMIP5 were used to climate change simulation. Since meteorological information of Mahabad synoptic station has a suitable time series, therefore, the statistics of this station were used for climate change analysis. Finally, we analyze the simulation results of climatic parameters for the future period.

Downscaling by change factor method
Using the delta method or the change factor, it is possible to calculate the desired climate variable changes in future periods compared to the previous period and the climate change scenario of that variable is obtained using the change factor method (Diaz-Nieto and Wilby 2005; Raju et al. 2018). In this case, the difference between the two variables is used for temperature. And for precipitation the ratio between them is used (Hulme and Jones 1994;Lipczynska-Kochany 2018). In the present study, these values are estimated for 2021-2045.
P∆ and T∆ represent the scenario of climate change, precipitation and temperature in the next period in each month, respectively. T future 20-year average temperature simulated by AOGCM in the next period for each month, T bases , is the 20-year average temperature simulated by the AOGCM in the observation period per month. P future and P base also represent the average of 20 years of simulated precipitation in the future and observed baseline, respectively. Thus, in the first step, the climatic variables of the main cell for the baseline and future were extracted from the IPCC site. Then, ArcGIS software was used to visualize the geographical coordinates of Cartesian. Then, by subtracting the long-term monthly values of the base period from the corresponding values in the next period, the climate change scenarios of each period were obtained. Finally, precipitation values were obtained by dividing the changes of future period values by the changes of base period values. Finally, the generated temperature and precipitation data under future scenarios are calculated from the following equation (Sunyer et al. 2012).
where T base and P base represent the time series of monthly observational temperature and precipitation in the base period, respectively. T and P represent the time series of temperature and precipitation resulting from climate change in the future and under a downscaling climate change scenario, respectively. In this study, climate change and scenarios of RCP26, RCP 45 and RCP 85 have been used to study the effect of climate change on the amount of inflow to the dam reservoir and its evaporation. To generate future precipitation and temperature data, three scenarios of the fifth IPCC report including RCP8.5, RCP4.5, and RCP2.6 were used. The specifications of the scenarios are given in Table 1. In each scenario, the effect of greenhouse gas emissions is classified based on its role on the level of inductions.

Performance evaluation criteria
Root-mean-square error (RMSE), mean absolute percentage error (MAPE), and coefficient of determination (R 2 ) were measured to assess the performance of scenarios and machine learning models used in this study. Mean absolute error (MAE) and correlation coefficient (CC) were also used to select the appropriate model for climate change (Arya Azar et al. 2022).
where x o is the observed (measured) value, x p is the predicted value, and n is the number of samples. x 0 and x p represent the average of the predicted and observed data, respectively.

Results of groundwater simulation
The groundwater level simulation was calibrated for steady and unsteady conditions. The trial-and-error method was utilized along with the PEST program in order to calibrate the model. In other words, first the range of changes in hydraulic conductivity and specific yield in each area was determined by trial and error; then, the PEST software was used to obtain the final calibrated values of hydraulic conductivity.
As depicted in Fig. 4, these values vary from 3 to 30 m/day. The lowest value of hydraulic conductivity, ranging from 3 to 3.5 m/day, belongs to the northern area of the aquifer. The mean hydraulic conductivity, ranging from 12 to 18 m/ day, was obtained in the center of the aquifer. The south of the aquifer is characterized by high hydraulic conductivity due to the presence of coarse-grained materials. Hydraulic conductivity in this area varies from 25 to 30 m/day. According to Fig. 4b, the calibrated specific yields range from 1 to 21% in unsteady conditions. The highest and lowest specific yields occurred in the southern and northern parts of the

c)
Observed Computed aquifer, respectively. The average range of this coefficient, varying from 12 to 16%, is related to the center and northeast of the aquifer. The trend of changes in simulated and actual groundwater balance under unsteady conditions in Fig. 5 reveals that these changes are different in different areas. For instance, in Fig. 5, the sinusoidal trend of changes shows that groundwater balance drops by about 1 m during the simulation period. Also, the time series indicates that MODFLOW has been able to correctly identify the trend of changes in groundwater balance in each month. Hence, it could be said that groundwater simulation has been properly carried out and the values of observation and simulation wells are very close to each other. Examining sudden changes in the actual and simulated values of groundwater balance in some steps (e.g., in Fig. 5-b and steps 11 and 12), we could infer that MOD-FLOW is capable of simulating sudden changes in groundwater and different stresses that are imposed on it.
Groundwater level was simulated to control the amount of groundwater withdrawal obtained from the optimization model results. The simulation was performed in two steady and unsteady conditions, and the results demonstrated an appropriate simulation performance for both conditions (Table 1). Further, to ensure the efficiency of groundwater simulation, we applied the model to the validation data, which resulted in RMSE, MAE, and R of 1.05, 0.90, and 0.98, respectively.

Integrated simulation of surface water and groundwater resources through the WEAP model
For the purpose of hydraulic simulation of the study area during the period 2006-2019, all the water supplies and uses (including the discharge of Mahabad River's headwaters, the catchment area of other rivers, agricultural, industrial, and urban consumption sectors, groundwater resources, and return flows) were separately entered as the input for the WEAP model. In order to minimize simulation error, we calibrated and validated the proposed model by using MARE Fig. 6 Simulated runoff and the actual values during the calibration and validation periods and NS indices on the basis of data from the hydrometric station of Mahabad. The first period (2006-2015) was considered for calibration, and the last 4 years (2016-2019) were used for validation of the model. The assessment results obtained from the WEAP model during the calibration period (NSE = 0.95, MARE = 0.41) and the validation period (NSE = 0.90, MARE = 0.61) for the hydrometer station indicate that this model has a great performance in the case of the station under investigation (Fig. 6). The trend of the simulated runoff time series as well as the values presented in Fig. 6 show that WEAP has successfully simulated the runoff. The level of accuracy of the calibration and validation periods is also acceptable.

Climate change results
The fifth report on climate change was used to study the effects of climate change on the rate of inflow and evaporation from dam reservoirs. The data of this report have different results in the production of precipitation and temperature data; however, selecting a suitable model that produces temperature and precipitation data in accordance with the data of the study area is of particular importance. In this study, 14 climate change models were used, and the best model for precipitation and temperature was selected using RMSE, MAE, and correlation coefficient. Thus, by comparing the historical data  of each model with the observational data, a model was selected using evaluation criteria. Among the available models, the EC-EARTH model, with RMSE, MAE, and correlation coefficient of 3.14, 1.8, and 0.99, respectively, had the highest agreement with the mean observational temperature data of the study area (Table 2). Therefore, this model was selected to extract future data (2026)(2027)(2028)(2029)(2030)(2031)(2032)(2033)(2034)(2035)(2036)(2037)(2038)(2039)(2040)(2041)(2042)(2043)(2044)(2045). The closest model to the EC-EARTH model is the GISS-E2-H model with RMSE, MAE, and correlation coefficients of 4.7.3, 5.5, and 0.78 respectively. Although other models such as CESM1 (CAM5) and IPSL-CM5A-MR had less error, their correlation values were about 0.83, which was less than the selected model. However, based on the minimum error values, the results of the FIO-ESM model were closer to the observational values. For this model, the RMSE, MAE, and correlation coefficients were 16.58, 15.4, and 0.4, respectively, which were lower than those of other models. For the precipitation parameter, the NorESM1-ME model with RMSE = 18.3, MAE = 15.1, and CC = 0.3 had the closest results to the selected model. After selecting the superior model, the trend of average temperature changes under three scenarios of RCP26, RCP45, and RCP85 was examined in the coming years. The results of temperature changes for the three scenarios are shown in Fig. 7. During 2026 and 2045, the average temperature will have a relative increase compared to the base temperature. According to the figure, most temperature changes will occur under the RCP85 scenario. For this scenario, the months of July to November will experience the highest temperature increase compared to similar months of the other two scenarios. In September, October, and November, the average temperature will rise by 2.45, 2.23, and 2.14 °C, respectively. The temperature in February will also rise from -0.71 to 0.39 °C. The lowest temperature increase will be under the RCP26 scenario. In the RCP26 scenario, there will not be much change in temperature for January and February, and it is observed that the temperature will decrease by 0.12 °C in February. Also, in the RCP26 scenario, the average temperature for April dropped by 0.22 °C. For the RCP45 scenario, the changes are not remarkable as RCP85 and are almost consistent with the RCP26 scenario. At RCP45, November will have the highest temperature rise (1.79 °C). The large temperature changes under the RCP85 scenario can be attributed to the increasing concentration of carbon dioxide. Moderate changes in average temperature can also be attributed to the constant concentration of carbon dioxide in the coming years. According to the RCP26 scenario, the concentration of carbon dioxide will generally decrease, which has led to small changes in the average temperature.
According to the selected model for precipitation (FIO-ESM), under the two scenarios of RCP26 and RCP45, in the first four months of the year and December of each year, the average precipitation will decrease. However, for the RCP85 scenario, precipitation will decrease almost throughout the year for the next period (2026)(2027)(2028)(2029)(2030)(2031)(2032)(2033)(2034)(2035)(2036)(2037)(2038)(2039)(2040)(2041)(2042)(2043)(2044)(2045). Nonetheless, in February, according to the two scenarios RCP26 and RCP85, precipitation will increase to some extent. In the RCP26 scenario, in May, June, and July, the precipitation will increase by 10.5, 6.92, and 2.7 mm, respectively, which among the studied scenarios, the most increase will belong to this scenario (Fig. 8). January and December, which are the rainy months of the study area, will experience a decrease in precipitation by about 10 mm under the RCP26 scenario. According to the RCP45 scenario, precipitation will decrease by 13.11 mm in January, which is the highest decrease in precipitation for a month among the three scenarios.

The impact of climate change on surface water and groundwater resources
Considering changes in precipitation and temperature as well as the balance of groundwater resources in the study area, we evaluated changes in surface runoff and discharge of the aquifer and applied them in the model. The results obtained from different climate change scenarios were incorporated to simulate the state of Mahabad aquifer using the WEAP-MODFLOW. Accordingly, the hydrograph of the aquifer was extracted of the climate change model under the average RCP 2.6, RCP 4.5, and RCP 8.5 emission scenario for the next 20-year period (Fig. 9). The results of climate change scenarios demonstrate that changes in groundwater balance have a negligible impact on groundwater resources. According to Fig. 9, the balance of groundwater varies from 0.8 to 15 cm in each month. The lowest reduction in groundwater balance occurred in January, while the highest reduction occurred in April, May, and June. In general, it can be concluded that the effect of climate change on the balance of wells in the region has been insignificant in almost all months.
The results establish that climate change will have little impact on groundwater balance in the future. Therefore, it can be concluded that climate change has positive impact on the reduction of groundwater balance and the amount of aquifer inflows. The findings of the present study also demonstrate that MODFLOW has strong performance in modeling the aquifer status and evaluating the impact of climate change scenarios on the amount of aquifer inflows. According to these results, groundwater resources of Mahabad currently have a good potential to be optimally and sustainably exploited in an integrated system along with surface water resources. The findings of Nematollahi & Sanayei (2022) are also consistent with the results of the present research. In addition, similar to the findings of Arya , our findings indicate that climate change in this area affects different parameters such as groundwater balance and precipitation (Soltani et al. 2023). In line with Dehghanipour et al. (2019), our assessments show that integrated simulation with WEAP-MODFLOW has strong performance in areas surrounding Lake Urmia. Therefore, this approach is suitable for integrated simulation of groundwater and surface water resources.
The simultaneous application of WEAP and MOD-FLOW models to Mahabad region confirmed that considering the integrated system of surface water and groundwater resources in areas that feature the interaction of an aquifer and a river facilitates simulating this interaction and estimating the volume of water exchange. The observed and simulated groundwater balance differed in various parts of the aquifer. There was a small difference between observation wells that were near the exchange point between the aquifer and the river. However, in the rest of the observation wells, there was a relatively large difference between the observed and simulated values, which implies the need for more time to calibrate this area. Despite the decrease in precipitation and increase in temperature, the impact of climate change on surface water and groundwater resources has not been substantial. Meanwhile, our evaluations showed that groundwater resources are more affected by climate change than surface water resources. These results are in line with the findings of Javadi et al. (2020) but in contrast to the report by Plunge et al. (2022), who found that climate change affects surface water resources more than groundwater resources. However, the continuation of the current trend and the possible effect of climate change on the water resources of the study area will cause a further decline in both groundwater and surface water. Therefore, future studies need to adopt appropriate management scenarios that could be implemented by the competent authorities to solve water problems. Unplanned water use in the agricultural sector, especially in dry areas where the only supply is groundwater, has led to a sharp decline in the water balance. In this context, presenting different solution scenarios without considering their effectiveness and evaluating the stability of the aquifer will not translate into better management of the aquifer. Furthermore, most studies have only focused on the effects of climate change and have not explored climate change adaptation, which needs to be addressed in future research.

Conclusion
Mahabad region, like many other regions, has aquifers with good potentials in terms of groundwater resources. However, it may be in decline in the future due to interactions between aquifers and rivers, overexploitation, and mismanagement of surface water and groundwater resources. Consequently, we used WEAP-MODFLOW in the present research in order to study the integrated performance of surface water and groundwater resources in Mahabad region and to examine the interactions between the two resources. We also aimed to study the impact of climate change scenarios of the fifth IPCC report on changes in groundwater balance. Therefore, we used five climate change models of the fifth IPCC report and three scenarios, namely RCP8.5, RCP4.5, and RCP2.6. The findings of this research demonstrate that the proposed model can implement an acceptable integrated simulation. The amount of return flows from the agricultural sector and the interactions between the aquifer and the river as obtained from the simulation were of acceptable accuracy and closer to the observational and actual values. The results of the integrated simulation of the study area show that interaction with the river has mostly influenced the southern and southwest parts of Mahabad aquifer, while exploitation mostly occurs in the northern part. Regarding the impact of climate change on the aquifer, the results demonstrate that the overall impact is not significant, while RCP 8.5 scenario shows more severe effects compared with the two other scenarios. Considering the impact of climate change on water resources, if the current trend continues and no effective changes take place, the water supplies of this region are likely to face a great challenge in near future. In order to avoid that situation, effective water management decisions must be taken. According to the findings of the present study, more research on this subject can help clarify other aspects of problems related to optimal exploitation of groundwater resources. Other researchers in this field can aim to study the effect of different exploitation strategies on the aquifer, the possibility of transferring water from Mahabad Dam to Lake Urmia, and the impact of changing cropping patterns in the region on groundwater resources.
Funding The author(s) received no specific funding for this work.

Conflict of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.