Climate Change Impact on Surface Water and Groundwater Recharge in Northern Thailand

Climate change is progressing and is now one of the most important global challenges for humanities. Water resources management is one of the key challenges to reduce disaster risk. In Northern Thailand, flood and drought have always occurred because of the climate change impact and non-systematic management in the conjunctive use of both sources of water. Therefore, this study aims to assess the climate change impact on surface water and groundwater of the Yom and Nan river basins, located in the upper part of Thailand. The surface water and groundwater regimes are generated by a fully coupled SWAT-MODFLOW model. The future climate scenarios are considered from the Representative Concentration Pathways (RCPs) 2.6 and 8.5, presented by the Coupled Model Intercomparison Project Phase 5 (CMIP5), in order to mainly focus on the minimum and maximum Green House Gas (GHG) emissions scenarios during the near future (2021–2045) periods. The results show that the average annual air temperature rises by approximately 0.5–0.6 ◦C and 0.9–1.0 ◦C under the minimum (RCP 2.6) and maximum (RCP 8.5) GHG emission scenarios, respectively. The annual rainfall, obtained from both scenarios, increased by the same range of 20–200 mm/year, on average. The summation of surface water (water yield) and groundwater recharge (water percolation) in the Yom river basin decreased by 443.98 and 316.77 million m3/year under the RCPs 2.6 and 8.5, respectively. While, in the Nan river basin, it is projected to increase by 355 million m3/year under RCP 2.6 but decrease by 20.79 million m3/year under RCP 8.5. These quantitative changes can directly impact water availability when evaluating the water demand for consumption, industry, and agriculture.


Introduction
Over the past few decades, the increase in global air temperature has caused several impacts, especially on climatic conditions. The number of extreme weather events has surged and there have been many natural disasters [1-3] because of the elevated CO 2 concentration and uncertainty in weather circulation, influencing changes in the hydrologic cycle [4][5][6][7]. As reported by the Intergovernmental Panel on Climate Change (IPCC) in the Fifth Assessment Report (AR5), the water cycle around the world has been changing and this change seems to be severe in the future due to the change of climate [8]. This causes an expectation that, in the future, flood and drought disasters will occur easily and severely in several regions [9][10][11][12][13], especially Southeast Asia, where adaptation of water management to avoid these natural disasters will become a major challenge in the current century [14]. Unfortunately, Thailand is one of the countries which is experiencing this issue. The intensity of floods and droughts increases every year and is expected to be more severe in the next 25 years, especially in the northern areas, for example, the Yom and Nan river basins. According to recorded history, the greatest flood disaster was in 2011 and caused damage to the northern and central areas [15]. Moreover, in 2015 and 2016, an extreme drought disaster occurred in the Yom and Nan river basins. This had an extensive effect on Thailand because it is an important agricultural and travel area, producing an enormous source of national revenue.
Referring to several pieces of research, one of the primary factors producing flood and drought phenomena is the amount of surface water and groundwater. Flood events can occur easily when the surface water cannot be drained immediately due to a large amount of groundwater in the wet season and limitation of capacity in the drainage system. In contrast, drought and water shortage can occur easily when there are less surface water and groundwater during the dry season. It can be seen that the change in climate conditions has a direct impact on the amount of surface water and groundwater, affecting the possibility of floods and droughts [16]. Therefore, the surface water and groundwater regimes must be simulated under different climate scenarios to prepare for various possible events. The widely used approach for indicating the projection of future climate is to apply a yield of Global Climate Models (GCMs) because of their reliable assumption in the evaluation of historical and future climate scenarios [17][18][19]. As an assumption of GCMs in the fifth phase of the Coupled Model Intercomparison Project (CMIP5), the future climate conditions are dependent on the Green House Gas emission scenarios, namely the Representative Concentration Pathway (RCP) [20,21].
Generally, the popular approach for studying the issue of climate change impact on the hydrological regime is the use of a hydrological model for simultaneous surface water and groundwater simulation. Nowadays, numerous hydrological models can be applied for the simulation of surface water and groundwater. SWAT and MODFLOW are two hydrological models, always linking with this issue [22][23][24]. However, the coupled usage of these two models is complicated due to their differences in spatial analysis, causing trouble for data preparation as well as model execution [25]. To solve this problem, the SWATMOD-Prep model was developed by Bailey et al., [4] with a graphical user interface (GUI). This platform consists of several functions for preparing and linking an input data file in order to spend less time and avoid errors during the settings process. Moreover, this coupled model is a semi-distributed hydrological model, considering only the important physical processes to provide an effective performance in the simulation of large scale watersheds over long periods of time. Prior to this study, the SWAT-MODFLOW model is applied to assess the conjunctive use of surface water and groundwater for irrigation of large-scale watershed, located in semi-arid regions [26]. It is found that the coupled SWAT-MODFLOW model provides satisfactory results regarding streamflow and groundwater elevations, which can be used for water resources management. In addition, Chunn et al. [27] used the SWAT-MODFLOW model to assess the climate change impact on the interactions of surface water and groundwater. The surface water and groundwater regimes are simulated individually. The generated surface water from the SWAT model is defined as an input in the MODFLOW model for a groundwater flow simulation. This approach can be called the semi-coupled simulation of SWAT-MODFLOW. In this study, the fully coupled SWAT-MODFLOW model is used to simulate the surface water and groundwater regimes simultaneously for a completed simulation of the hydrologic process of the watershed. Furthermore, the effective replacement of the groundwater module in SWAT by MODFLOW can provide the output of surface water and groundwater simulation at a finer resolution than using only the original SWAT model or semi-coupled SWAT-MODFLOW model. Therefore, this study aims to assess the climate change impact on surface water and groundwater recharge (percolation) of the Yom and Nan river basins, located in the north of Thailand. The trends of both sources of water under each future scenario can be applied to anticipate the intensity of drought disasters that may occur in the future and could also be an approach for the development of the water management plan in the near future period.

Hydrological Model
The hydrological model mainly used in this study is a fully coupled model of SWAT-MODFLOW. The original SWAT groundwater part is replaced by the groundwater flow process of the MODFLOW model, simulating alongside with the land surface section of the SWAT model as one system. The Hydrologic Response Units (HRUs) and sub-basin-based variables of SWAT are linked to grid-based variables of MODFLOW. The recharge rate from SWAT's HRU is mapped to the Recharge package of MODFLOW; then, the stream-aquifer interaction is simulated and linked back to the sub-basin channel of SWAT. The HRUs are disaggregated and intersected with MODFLOW's grid to estimate an interaction between surface water and groundwater flows. Meanwhile, MODFLOW's river cell is intersected with SWAT's sub-basins to indicate the interchange of the river stream and aquifer. During the simulation, the groundwater process of the SWAT model is replaced by the MODFLOW module, performed as a subroutine of SWAT in the temporal resolution of a daily time scale [28]. The linkage parameter between both models is the water percolation of each subbasin (SWAT) and river recharge from the aquifer (MODFLOW). The land surface flow section of the SWAT model depends on the water balance equation following Equation (1), while the streamflow section is represented by Manning's equation (Equation (2)) for estimating a river discharge as free surface flow in an open channel and routed by the Muskingum River routing method as Equation (3) [29]. The groundwater regime in the MODFLOW model is simulated as a spatial discretization with grid blocks. The movement of groundwater is generated in the finite-difference form following Equation (4) [30].
q out,2 = C 1 · q in,2 + C 2 · q in,1 + C 3 · q out,1 , where SW 0 is the water content of soil on the first day (mm.), SW t is water content of soil on the last day (mm.), P is the amount of precipitation on day i (mm.), Q Land is the surface runoff on day i (mm.), Q River is river runoff on day i (mm.), E is the amount of evapotranspiration on day i (mm.), W is percolation exiting the lower soil profile on day i (mm.), Q gw is the amount of water recharge on day i (mm.), A is the channel cross section (m 2 ), R is the hydraulic radius (m.), S is the channel slope length, n is the Manning's coefficient of river channel (m.), q in,1 is the inflow rate at the first-time step (m 3 /s), q in,2 is the inflow rate at the last time step (m 3 /s), q out,1 is the outflow rate at the first time step (m 3 /s), q out,2 is the outflow rate at last time step (m 3 /s), C is the Muskingum coefficients, Q i is the groundwater flow rate (m 3 /s), SS is Specific Storage (m -1 ), ∆V is grid cell volume (m 3 ), and ∆h is change of water head over a time interval of ∆t.

Study Area
The Yom and Nan river basins are important and vulnerable regions located in Northern Thailand between 14 • 50 N and 99 • 16 E, and 18 • 37 N and 101 • 21 E (Figure 1). This area covers approximately 58,782.93 km 2 or 37% of the Chao Phraya river basin. The major land use in this study is agriculture, especially paddy and plant farms. Consequently, water availability is mainly used for agriculture. The main rivers are the Yom and Nan rivers flowing along the watershed from the steep valley in the upstream area through the narrow plain area and terraced mountain in the middle until the flat plain area downstream. The flood and drought disasters have occurred in both basins several times due to a different river capacity between upstream and downstream areas, lack of rainfall, and unsystematic water management, even though there is the Sirikit Dam in the Nan river basin. With regard to hydro-climatic information, Figures 2 and 3 provide the information concerning an average monthly runoff, rainfall, maximum and minimum temperatures of the Yom and Nan river basins, respectively. This information is an average value, recorded over a period of 31 years from 1985 to 2016. It is noticed that the Yom and Nan river basins have a similar pattern in this hydro-climatic information. The large majority of rainfall occurs during the wet season (May-October). This is due to the effect of depression, the Southwest monsoon, and the Northwest monsoon [31]. The highest rainfall intensity is in September at 250 mm, approximately, while the lowest is in January at 3 mm. The peak runoff happens in October with around 23,000 and 45,300 million m 3 for Yom and Nan river basins, respectively. The average maximum air temperature increases during January and hits a peak in April at around 31.3 • C, then it is reduced until December at 25 • C. For the average minimum air temperature, it rises from January (18 • C) to April (25 • C), as same as the average maximum air temperature, and slightly decreases to 23.5 • C at October before sharply dropping in November (21 • C) and December (17 • C).

Data Collection
The required data used for implementing can be classified into three main types. First, the geological data represent the basic physical characteristic of the watershed, consisting of the topography, land use, and soil property. The surface elevation was defined by a 90 m resolution Digital Elevation Model (DEM) produced by the Shuttle Radar Topography Mission (SRTM) of the National Aeronautics and Space Administration (NASA). Aquifer thickness was defined by the soil depth data from the global gridded soil information, produced by the International Soil Reference and Information Centre (ISRIC) with a resolution of 1 km. The land use and soil properties data were provided by the Land Development Department (LDD) and the Department of Mineral Resources (DMR), respectively. The land use data, created in 2016, is divided into five main types for representing different growing seasons which are the forest, crop field, fruit farm, plant farm, urban and water resource. The soil properties data are collected in 2016 and categorized into 15 soil types based on Thai soil types by DMR.
Second, the hydrological data were a time-series required for model calibration and definition of water operation of the dam. These data consisted of the river discharge, groundwater level, and reservoir operation. The river discharge data is recorded daily from a runoff station operated by the Royal Irrigation Department (RID). The reservoir operation data is a daily outflow release of the Sirikit Dam, provided by the Electricity Generating Authority of Thailand (EGAT). The groundwater elevation is monthly-measured data from observation wells operated by the Department of Groundwater Resources (DGR). The location of runoff stations, observation wells, and the Sirikit dam are displayed in Figure 1.
The third type of data was meteorological, indicating the weather conditions of the watershed such as air temperature and rainfall. These data are separated into two parts for individual simulation, consisting of the Global Climate Model (GCM) data and observation data. The GCMs data are output from the Coupled Model Intercomparison Project Phase 5 (CMIP5), which takes weather conditions generated on a global scale and divides them into historical and prediction periods. These data were provided by the National Institute for Environmental Studies, Japan (NIES). The names of GCMs used for future scenario simulations in this study are listed in Table 1. The observed climate data were collected from the rainfall and weather stations of the Thai Meteorological Department (TMD) in the locations shown in Figure 1. Only rainfall data is used on a daily scale, while the other variables such as the air temperature, wind speed, and solar radiation are considered on a monthly scale because of the data available.

Model Set-Up and Evaluation
The surface water and groundwater regimes were generated via SWATMOD-Prep, which is a preparation interface for linking SWAT and MODFLOW models. Initially, the simulation is acted through the ArcSWAT2012 model in a temporal resolution of daily, over the period of 2004-2016 by defining 2004-2006 as a warm-up period. The watershed is divided into 111 sub-basins by depending on DEM and covering most areas of the Yom and Nan river basins. There are two outlet locations, defined at station Y.5 and N.8A for the Yom and Nan river basins, respectively. The HRUs are classified into 2102 numbers by the overlay of sub-basin, land use type and soil type in a definition ratio of 1% for these three variables. The reservoir operation of the Sirikit dam is indicated by using daily data of measured water release along the simulation period. Then, the constructed SWAT model files are used to create and link with the MODFLOW files by the SWATMOD-Prep, in which the spatial resolution is 4 km due to a limitation of grid cell numbers of this version. The aquifer parameter is defined as a raster, based on the soil types. The period and temporal resolution of the SWAT-MODFLOW simulation are followed, setting in the SWAT model. With regard to an output, the water yield, generated by the SWAT-MODFLOW model, is a combination of surface runoff (SWAT), lateral flow (SWAT) and river recharge from the aquifer (MODFLOW). So, in this study, the water yield would be considered as surface water in order to estimate the water availability of the ground surface. Meanwhile, the groundwater recharge can be estimated by the water percolation of surface flow from SWAT only because the other Recharge package of MODFLOW is not included in this study.
The constructed model was investigated in the calibration process to find the most appropriate value for each parameter. The simulated results, such as stream discharge and groundwater elevation, were compared with the observation data from 2007-2011 for calibration and 2012-2016 for validation. In addition, the model efficiency coefficients, namely, the Coefficient of Determination (R 2 ), Nash-Sutcliffe Efficiency (NSE) [32], Root Mean Square Error (RMSE) and Percent Bias (PBIAS) [33] were applied for model performance investigation. According to several kinds of research, these coefficients are widely used and can be computed following Equations (5)- (8). The value of NSE is between -∞ and 1, while R 2 is between 0 and 1. The model will provide high performance when the NSE and R 2 values are closer to 1. The RMSE and PBIAS values will be closer to 0 if the model has high accuracy because this means that there is a small difference between the observed data and simulated result. The positive and negative values of PBIAS mean that the simulated results are under-or overestimation, respectively [33,34].
where Y i and Y ' i are the observation data and simulated result, respectively, and Y is the observation data on average. n is the total number of considered data.

Bias Correction Method
Initially, the GCMs data were evaluated by comparing them with the observation data during the historical period in order to select an appropriate GCM, providing a consistent value and pattern with the observation data. This is because each GCM produces a different output due to the diversity of assumptions and methods applied. However, there is still a systematic distributional bias because the GCMs data are generated on a global scale that differs from the resolution of the observational data. Therefore, in this study, the selected GCMs data had to remove this bias by using the statistical technique, namely, bias correction methods. It is generally known that there are several ways to implement bias correction methods, but the most widely used option is the Change Factor Method (CFM) because of its simplicity and straightforwardness. The approach of CFM can be applied by basing it on the characteristic of climatological data. The average monthly, seasonal or annual values are better than the daily value because a higher frequency of temporal scales probably causes a decrease in a GCM's reliability [35]. The Change Factors (CFs) [36] were computed in this study for each month of the year separately from a difference or ratio between the original GCMs data during the historical and projection periods. These CFs are added or multiplied to the daily observation data, following Equations (9) and (10) for air temperature and precipitation data, respectively [37,38]. This can be called the Shifting and Scaling method [39]: P bias−cor. y,m,d where T is air temperature ( • C) and P is precipitation (mm.). The upper bar is the average value for the period. Superscripts bias-cor., obs, and org are bias-corrected, observed, and original GCMs values, respectively. Subscripts future, baseline, y, m, and d are prediction and reference values in the annual, monthly and daily time scales, respectively.

Assessing the Impact of Climate Change
The projected historical observation corrected based on CFs were inputted into a calibrated SWAT-MODFLOW model for the future scenarios simulation of surface water and groundwater. For assessing the surface water and groundwater changes, the simulated results of the reference period and projected period (2021-2030) were compared. The future scenarios depended on the global annual Green House Gas (GHG) emission scenarios proposed in the Coupled Model Intercomparison Project Phase 5 (CMIP5), as a Representative Concentration Pathway (RCP) of 2.6 and 8.5 to focus on the minimum and maximum conditions of GHG emission issues, respectively. The climate change impact on surface water and groundwater was analyzed by comparing the air temperatures and rainfall with variations in both water resources. Moreover, the amount of surface water and groundwater from each scenario was correlated with an expected water demand to assess water availability under divergent climate scenarios.

Model Evaluation
The comparison of surface water (river runoff) and groundwater (groundwater elevation) from the simulated results and observation data are presented in Table 1 and Figures 4-7. Table 1 provides the information about the NSE, R 2 and RMSE values used to evaluate the river discharge and groundwater elevation, obtained from simulation and observation. It can be seen that during the calibration period (2007-2011) the NSE and R 2 values of the runoff station are greater than 0.75. The performance of the model simulation seems to be satisfactory because of the NSE, R 2 value close to one when considered on a daily basis. The NSE with the nearest value to 1 for the Yom and Nan river basins is at stations Y. 16 and N.60, with values of 0.857 and 0.842, respectively. Similarly, the highest R 2 values for the Yom and Nan river basins were 0.878 and 0.874 at stations Y.16 and N.60, respectively. For the investigation of groundwater, it is necessary to evaluate a results from all observation wells holistically and displayed as quartile plots because the groundwater elevations data is observed and recorded on a monthly basis, so there are too few observation data to evaluate using the efficiency coefficient index in the individual observation wells. The results show that the NSE and R 2 values of all observations wells are 0.710 and 0.765, respectively. Moving on to the results of the validation period (2012-2016), most NSE and R 2 values of runoff stations are lower than received from the calibration period. This might be because the calibration period was included by many flooded years, having a large amount of river discharge, while almost every year in the validation period has a low discharge value. Therefore, the validated results are probably overestimated when using the same parameters from the calibration period. However, all of them (NSE and R 2 ) are of satisfactory value. Whereas, the NSE and R 2 values of groundwater elevation during the validation period are 0.761 and 0.829, respectively, which are higher than received from the calibration period. In addition, the RMSE values are a possible minimum value for each location that can be simulated by the SWAT-MODFLOW model. The PBIAS values of streamflow are between -10.081% and 10.351% in the calibration period, and between −35.597% and 14.293% in the validation period. It can be noticed that most of the PBIAS values obtained from both calibration and validation periods are negative values, meaning that the streamflow from the SWAT-MODFLOW model is an overestimation. The PBIAS values of the stations Y.16 and Y.5 show a large negative value in the validation period of approximately −35%. This is because the amount of water in the validation period is very low. Thus, the summation of discharge values in the validation period is small and causes a large value of the PBIAS index when using to divide.
While, for groundwater, the PBIAS values of groundwater elevations are defined as −10.630% and 12.363% in the calibration and validation periods, respectively. This means that the groundwater elevations are overestimated in the calibration period and underestimated in the validation period.
The temporal comparison between the observed data and simulated results of surface water and groundwater are shown in Figures 4-7. Figures 4 and 5 present the comparison of daily discharge from observation data and simulated results at runoff station Y.16 and inflow of the Sirikit dam. Figure 6 presents the comparison of monthly groundwater elevations from observation data and simulated results at observation well NT.57. These figures display that the values from SWAT-MODLOW simulation are consistent with observation data during both calibration and validation periods. With regard to Figure 7, showing a comparison of groundwater depth during the calibration (a) and validation (b) periods, the groundwater elevations from all observation wells were converted to the groundwater depth and plotted as a quartile plot in order to avoid an excessive differentiation in the elevation of observation wells, which can cause an unclear result for the reader. It proves that this model can provide the pattern of groundwater levels similar to those recorded from the observation wells. This is because the large majority of results showing in the 1:1 comparison graph is close to the perfect agreement line. This indicates that the groundwater level from simulation data is consistent with the observation data.                   Table 2 represents the quantitative consistency of monthly rainfall intensity between the original GCMs data and observation data through the NSE and RMSE values. The monthly rainfall intensity is averaged from the data between 1980 and 2005 due to the end of the historical period of GCMs. Three GCMs provide the NSE values in good criteria (higher than 0.8), namely, MIROC5, MPI-ESM-MR, and CNRM-CM5. The MIROC5 provides the greatest value at 0.93, which is slightly higher than CNRM-CM5 (0.9) and MPI-ESM-ES (0.86). Besides, these three GCMs also show lower RMSE values than the others at 19. HadGEM2-ES, which was underestimated. Therefore, CNRM-CM5, MIROC5, and MPI-ESM-MR were selected for use in the future scenarios simulation of surface water and groundwater regimes in Yom and Nan river basins.

Climate Change Scenarios in Yom and Nan River Basins
The analyzed results of climatic conditions in Yom and Nan river basins are presented in  Figure 9b, the annual rainfall of the Yom and Nan river basins has a consistent trend with the air temperatures and there is a possibility that it will rise in the near future. This is because the average value of annual rainfall from most future simulation cases is higher than the reference period under both minimum (RCP 2.6) and maximum (RCP 8.5) GHG emission scenarios. The average annual rainfall in the Yom river basin grows from around 1140 mm at the reference period to 1160-1240 mm under scenario RCP 2.6 and 1190-1120 mm under scenario RCP 8.5. The average annual rainfall in the Nan river basin grows from around 1260 mm in the reference period to 1340-1430 mm and 1340-1400 mm under scenarios RCP 2.6 and RCP 8.5, respectively. However, it is not a completely upward trend because there is an exception from the analyzed result of MIROC5 in the Nan river basin, showing a reduction of average annual rainfall from the reference period at 1260 mm to 1220 and 1240 mm during the near future under scenarios RCP 2.6 and RCP 8.5, respectively. Figures 10-12 demonstrate the change of rainfall and air temperature in the monthly variation. It can be seen that both maximum and minimum air temperature of the Yom and Nan river basins are raised in the future under scenarios RCP 2.6 and 8.5. Under scenario RCP 2.6, there are two periods of the year, showing an increase in the maximum air temperature which are March-May and September-November. The maximum air temperature can increase by 1-1.5 • C from the reference period at those times. While, for the scenario RCP 8.5, the maximum air temperature is totally increased all year round. The highest increase mostly occurs in October where the temperature can rise by 2 • C from the reference period. Similarly, the minimum air temperature has the possibility of increasing throughout the year in the future under both RCP 2.6 and 8.5, especially during April-September. With regard to the change of monthly rainfall intensity, the Yom and Nan river basins show a different trend for this matter. In the Yom river basin, under scenario RCP 2.6, the rainfall seems to be increased from the reference period during May-July around 10 mm/month (on average) and decreased at the same rate in September. This is almost similar to what happened under RCP 8.5 where there is a 20 mm increase in monthly rainfall in June, but it equals to the reference period in the other month of the year, not decreasing in September, which is the same as under RCP 2.6. In the Nan river basin, the results of CNRM-CM5 and MPI-ESM-MR show that the amount of monthly rainfall intensity is increased (approximately 25 mm/month) during May-September in the overview under both scenarios RCP 2.6 and 8.5. Whereas, the results of MIROC5 are that the amount of rainfall is decreases in May for 40 and 60 mm under RCP 2.6 and 8.5, respectively, but other months are similar to the reference period.

Impact of Climate Change on Water Availability
The annual water yield and percolation during 1981-2005 (reference period) and 2021-2045 (near future period) are illustrated in Figures 13-17 by comparing to the amount of water demand, reported by the Department of Water Resources. According to Figure 13, based on the water demand of the Yom and Nan river basins, amounting to 3575 and 4638 million m 3 per year, respectively, it can be considered that there are 11 water shortage years (red bar) and two years that needed percolated water for the water supply (orange bar) in the Yom river basin. In the Nan river basin, there is neither a water shortage year nor a year that needed percolated water. This might be an influence of water management operated by the Sirikit dam. The most lethal water shortage year obtained from the simulation is 1993, which corresponds to what has been recorded in Thailand. Moreover, in this year, the impact of water shortage had severely and widely damaged the Northern and Central regions of Thailand. Regarding the future scenarios, the result from all three GCMs displayed in Figures 14-17 shows that the amount of surface water (water yield) in the Yom river basin decrease to 393.85 and 282.71 million m 3 under the minimum (RCP 2.6) and maximum (RCP 8.5) GHG scenarios, respectively. While, in the Nan river basin the surface water projected under scenarios RCP 2.6 and 8.5 is increased 347.67 million m 3 and decreased 12.12 million m 3 , respectively. Similarly, the quantity of groundwater recharge (water percolation) in the Yom river basin under scenarios RCP 2.6 and 8.5 is expected to decrease by 50.13 and 34.07 million m 3 , respectively. There is a 7.34 million m 3 increase and an 8.97 million m 3 decrease in the groundwater recharge of the Nan river basin under scenarios RCP 2.6 and 8.5, respectively. Besides, the number of water shortage years in the Yom river basin will increase significantly under both the minimum (RCP 2.6) and maximum (RCP 8.5) GHG emission scenarios. The number of years with water percolation required will increase under scenario RCP 2.6 but decrease under scenario RCP 8.5. In the case of CNRM-CM5 and MPI-ESM-MR, the number of water shortage and percolated water required years are one above the reference period under scenario RCP 2.6. Meanwhile, MIROC5 shows that the water shortage years increase from four years from the reference period, but the percolated water required years remain constant. However, under scenario RCP 8.5, these three GCMs indicated that the number of years with water percolation needed dropping one year from the reference period, while the number of water shortage years increased by 2-3 years. In addition, there are at least two years in the near future period with the possibility of a more severe water shortage year than in 1993. The number of water shortage and percolated water years needed in the Nan river basin is still zero even though there are several years where the amount of water yield and percolation are decreased to around the same level as the amount of water demand.   1982  1983  1984  1985  1986  1987  1988  1989  1990  1991  1992  1993  1994  1995  1996  1997  1998  1999  2000  2001  2002  2003  2004 1982  1983  1984  1985  1986  1987  1988  1989  1990  1991  1992  1993  1994  1995  1996  1997  1998  1999  2000  2001  2002  2003  2004

Discussion
According to the model efficiency coefficients estimated from the calibration method, it is clear that the SWAT-MODFLOW model provides satisfactory performance in the simulation of surface water and groundwater regimes. The NSE and R 2 received from the comparison of simulated and observed river discharge data on a daily scale are higher than 0.75, which is defined as good criteria. This is consistent with findings in [40], that is, the SWAT model shows congenial results in a simulation of surface water regimes in the Chao Phraya Watershed (Yom and Nan river basins are included). The value of NSE and R 2 in their simulation is greater than 0.5. For the groundwater regime simulated by the MODFLOW module, the comparison of groundwater depth from the simulation and observation in the 1:1 graph has shown that they are closed, which is similar to what was found in [41]. The groundwater level from their groundwater simulation using the MODFLOW model was close to the observation data, even though the calibration period is different from this study. Moreover, based on the comparison of rainfall intensity from the original GCMs output and observation data, this study found that the most appropriate GCMs used to project the future climate scenarios in the Yom and Nan river basins are MIROC5, CNRM-CM5, and MPI-ESM-MR. This is consistent with several studies that investigated climate change in Thailand and used the GCMs output from CMIP5 to define future climate conditions. For example, Kamworapan and Surussavadee [42] evaluated the performance of GCMs from CMIP5 for the projection of climatological temperature and precipitation. The historical simulation of climatological conditions from 40 GCMs were compared in their study and they pointed out that these three GCMs are in a good classification, producing a negligible error when compared with observational data. Notably, CNRM-CM5 is one of the six GCMs that provide the best performance for both temperature and precipitation.
Pratoomchai et al. [43] assessed groundwater availability in the Upper Chao Phraya river basin (Yom and Nan river basins are included in this area) under climate change impacts throughout 2026-2040. Their study found that the annual rainfall of MIROC increases 0.1-6.1% from the reference period (1986-2000) on average. In addition, there was a significant upward trend in the air temperature, approximately 1.3-1.8 • C. The maximum GHG emission scenario (RCP 8.5) provides the highest increase in air temperature between 2.5 and 3.8 • C. Meanwhile, in this study, an increase of rainfall, minimum and maximum air temperature was found during 2021-2040 in both the Yom and Nan river basins. However, the average annual air temperature increased 0.5-0.6 • C under the minimum (RCP 2.6) GHG emission scenario and 0.9-1.0 • C under the maximum (RCP 8.5) GHG emission scenario, which is lower than that estimated in the study of Pratoomchai et al. [43]. The highest increase in the average annual minimum and maximum air temperatures is approximately 1.34 and 1.02 • C, respectively, as shown in the simulation case of MPI-ESM-MR under scenario RCP 8.5. The annual rainfall increased by approximately 20-200 mm, or 1.6-16%, from the reference period , which is higher than the result found by Pratoomchai et al. [43].
The assessment of surface water (water yield) and groundwater recharge (percolation) under various future climate scenarios implemented in this study is at a regional scale. Consequently, the bias correction method has to apply with the future climate condition from GCMs in order to remove a gap between GCMs output and observational weather data. The simulation of future scenarios directly using the single GCM runs as input is quite a preliminary approach. This would be a source of uncertainty influencing the climate condition evaluated. The intensity of air temperature and precipitation might be underestimated when compared with a projection on a global scale. Usually, hydrologists make use of RCM outputs and even these are sometimes too coarse for hydrologic predictions. However, in this study, future simulations of surface water and groundwater regimes are implemented during the near future period (2021-2040) in order to avoid the effect of land use change. Hence, these analyzed results can be used as the approach for future study related to this area, and also in planning for an elementary water management plan. Based on this consideration, if the amount of GHG emission is at a high rate continuously and tends to be scenario RCP 8.5, it is necessary to plan for a more scrupulous policy of water usage. This is because there is a high risk that the water availability will be decreased from the ground surface or underground, due to the increase in air temperature and a lower amount of precipitation.

Conclusions
As the calibrated results of the R 2 , NSE, RMSE, and PBIAS indexes are classified as good criteria in several locations, especially in the calibration period, the coupled SWAT-MODFLOW model has proved that it can be applied to generate the surface water and groundwater regimes in the Yom and Nan river basins satisfactorily. Besides, there is an obvious upward trend in rainfall, and minimum and maximum air temperature from the reference period to the near future period in both the Yom and Nan river basins. The average annual air temperature increases by approximately 0.5-0.6 • C under the minimum (RCP 2.6) GHG emission scenario and 0.9-1.0 • C under the maximum (RCP 8.5) GHG emission scenario. Meanwhile, the annual rainfall increases by approximately 20-200 mm, but there are no differences between the minimum (RCP 2.6) and maximum (RCP 8.5) GHG emission scenarios. Therefore, based on the similar increasing trend of the annual rainfall under both scenarios, the amount of water yield and percolation under scenario RCP 8.5 became lower than that received from the reference period and scenario RCP 2.6 because of the higher air temperature. In the Yom river basin, the years with percolated water needed are almost constant compared to the reference period under scenario RCP 2.6. While under the scenario RCP 8.5, it was found to be a water shortage year due to a greater lack of water availability. However, the water shortage year significantly increased during the near future under both GHG emission scenarios. Moreover, there is a possibility that the water shortage year during the near future period can become more severe than any that occurred in the reference period. Whereas, the water shortage and water percolation needed years of the Nan river basin remains constant from the reference period. This means that the Yom river basin is expected to face more drought events under both scenarios RCP 2.6 and RCP 8.5. There is no concern for the Nan river basin due to the benefit coming from the Sirikit Dam. However, the intensity of drought events is more severe if the GHG emission remains at a high level or in the maximum scenario (RCP 8.5).

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