Water Scarcity Risk Assessment for Multi-Administrative Units in Agricultural Watersheds Using Integrated QSWAT-WEAP and GIS-Based Approach

1. Introduction

Climate change has significantly affected hydrologic cycles and water resource management over the past decade. Particularly, an increase in global average temperature leads to the changes in rainfall, runoff, and increased evaporation rate [1]. These changes cause higher severity of water-related disasters especially unusual prolonged drought in many regions around the world during 2020-2023 [2]. Impact assessment and a plan for coping with the climate change require a numeric model with high precision to simulate hydrological processes and prediction of global changes [3,4]. Southeast Asia is one region that has been most significantly affected by the climate change. Cambodia, Vietnam, and Laos, in particular, has been affected by the drought during 2020-2022, impacting region’s agricultural productivity and food security [5]. Water management in this region is complicated since both physical and socioeconomic factors have to be considered. As a result, models for hydrologic and water management become more attracted [6], and this results in the development of decision support system integrating several models for systematic problem analysis and strategic planning in water management [3,7].

Thailand has experienced severe drought over the past decade, especially the Northeast where agriculture is damaged over 15,000 million Baht in 2022 [8,9]. Applying hydrologic models to water prediction and management planning is unlimited in Thailand [10]. Most of the previous studies have tried to use single models for problem analysis, such as using hydrologic models to assess the impacts of land use changes on the runoff in the Northeast [11], using water budget models to analyze water demand in the Chao Phraya River Basin [12]. However, studies on integration of several models are not sufficient although the potential of problem analysis is high [7]. Moreover, those studies emphasize watershed or province level analysis, so the analysis at smaller level for supporting local planning is still lacking [13]. Development of spatial water assessment and monitoring system is important to improving the efficiency of country’s water resource management [6], especially a system to support decision making at local level integrating with national water management system [7,14]. Integrating hydrologic model technology with geographic information system is essential for decision support system development on future water management [15].

Currently, spatial model is employed in situational analysis of water because it can display geospatial and temporal data to support stakeholders to have more understanding of situations easily [16]. Integration between hydrologic model and water budget model is more widely used in many regions throughout the world for spatial water prediction. In Cambodia, for example, Touch et al. [17] used both of these models to assess the impacts of urban development on water security in the Tonle Sap Lake Tributary Basin. The study of Tran et al. [16] shows the quality of water management planning for semi-arid areas in Cuba. Also, in Vietnam, Khoi et al. [18] applied both mentioned model to assess water scarcity risk at sub-basin level in Vietnam.

One of the hydrologic models globally used for hydrologic analysis in river basins is SWAT because it can be used to simulate complicated hydrologic processes at sub-watershed level [19]. A study conducted by Bal et al. [20] demonstrates that SWAT can be used to simulate the impacts of land use changes in India on runoff precisely. This is consistent with a study conducted in Turkey by Aibaidula et al. [21] using SWAT to analyze the effects of climate change on reservoir water storage volume. In Ethiopia, Eshete et al. [22] used SWAT to analyze water scarcity risk in drought-prone areas. The important advancement of SWAT is that it was developed into a public domain model that can help researchers develop their model without financial constraints. This can lead to the integration of Quantitative Soil and Water Assessment Tool (QSWAT) [23,24,25] and Quantum Geographic Information System (QGIS) [25,26] to increase capacity in spatial analysis. Currently, QSWAT has been widely used for hydrologic model development in Southeast Asia, such as an analysis related to the effects of climate change on rained-dependent agricultural sites in Asia conducted by Goswami et al. [27], a study on geospatial and temporal hydrologic elements in a basin at regional level in Thailand [28,29,30], and potential assessment of water resources for dry-season crop cultivation planning [31].

Also, Water Evaluation and Planning System (WEAP) is a model widely used because it has the strengths in analyzing water balance and water scarcity in a drought-prone area, and it has been widely applied in many regions both in historic-present and future [32,33]. In South Africa, WEAP was utilized by Al-Mukhtar and Mutar [34] to analyze water allocation for different sectors in dry season while Ben Salem et al. [35] applied it to irrigation system and water scarcity risk assessment in drought-prone areas in Morocco. In the context of Southeast Asia, WEAP is also commonly used. In Vietnam, for example, the research conducted by Dau et al. [36] applied it to the water scarcity risk assessment for irrigated area in the Mekong Delta. In Thailand and Southeast Asia, Kuntiyawichai and Wongsasri [37] assessed the impacts of drought on water security at community level while Forni et al. [38] and Guzman et al. [39] developed a decision support system for water management in agricultural sites with water scarcity risk in Cambodia and Philippines, focusing on community participation in water utilization planning.

However, the previous studies lack small-scale spatial analysis, especially at local government organization level [40,41]. The study conducted by Candido et al. [42] and Nagata et al. [43] report the challenges of the development of output system that is easily understood by local users. Thailand has encountered several constrains in local water management, such as a lack of efficient tools for analyzing and forecasting water situation and limited data sharing between central and local government organizations, which cause inefficiency in local drought management. Therefore, this research emphasizes integration of QSWAT and WEAP models to solve these mentioned constrains through a developed sub-district level water scarcity risk assessment that can display the output in clear maps and graphs. Moreover, an open-source model was used for system modification to fit the local context. Integration of both models also increase capacity in hydrologic analysis and water resource management [33], which are essential for water management planning in drought-prone areas in Thailand.

This study aims to integrate QSWAT and WEAP models for sub-district level water scarcity risk assessment in Maha Sarakham Province, which is regarded as a drought-prone area in the Chi River watershed, northeastern Thailand. Maha Sarakham has repetitively faced severe drought over the past decade, especially in a prolonged dry season that agricultural sector is affected. Although the Chi River is the major water resource in this province, allocation of water into remote areas is limited because agricultural sites and communities have been expanded. As a result, QSWAT-WEAP integration is the key to systematic water situation analysis in a drought-prone area particularly in considering streamflow at sub-watershed level using QSWAT together with assessing sub-district level water scarcity risk using WEAP. Furthermore, this study focuses on developing decision support system that can display clear maps of sub-district level water scarcity risk in order that local government organizations and stakeholders can use it as a tool for planning their drought solutions in drought-prone areas, where water demand is high but water source is limited.

The result of this study can be useful for any agencies responsible for water management and drought solution at both policy and operational levels, such as Office of the National Water Resources, Royal Irrigation Department, and Local Government Organization, in order to assess water scarcity risk and establish mitigation measures for drought-prone areas. The displayed output of low risk can allow local administrators and officers to understand situations and provide high efficiency of resource allocation policy especially in Maha Sarakham, where water management for drought-prone area is very challenging. Integrating QSWAT and WEAP models with Geographic Information System (GIS) can improve the precision of analysis and forecast of water scarcity risk. Farmers and water consumers can utilize data from this system for their cultivation planning and water consumption aligned with water source available at each time period. In addition, the methodology of the study can serve as a model for system development of water scarcity risk assessment in other drought-prone areas, where water allocation is limited and water demand is various.

2. Materials and Methods

2.1. Study Area

2.1.1. Hydrological Characteristics of the Study Area

The study area is within the Chi-Mun Rivers system in northeastern Thailand (see Figure 1) including diverse topographic features but mostly the Korat Plateau with elevation from 193 to 371 meters above mean sea level (MSL). On the western side, the topographic features include mountains in Phetchabun Province with elevation from 729 to over 1,622 meters above MSL. Its hydrologic network is composed of two main rivers including the Chi River, where its source is from a western mountainous area, and the Mun River, which flows through the southern part of area. Additionally, it consists of two large reservoirs; Ubolrat and Lam Pao Reservoirs, which help control streamflow in the watershed. There is also a network of six meteorological monitoring stations (E.22B, E.91, E.9, E.66 and M.95) strategically distributed throughout the watershed to monitor spatial variability of flow patterns. Diverse physical features and comprehensive inspection network of this area become significantly reasonable to develop hydrologic modeling by using QSWAT.

2.1.2. Administrative Boundaries and Adjacent Areas

The major study area is Maha Sarakham Province including 5,292 square kilometers (sq.km) located in a strategic location, where several tributaries within the Chi-Mun system joining together. The Province borders Khon Kaen Province to the north, Roi Et Province to the east, and Surin Province to the south. It consists of 133 sub-districts, and each of them is clearly defined boundaries for systematic assessment and water resource management. Nevertheless, this study area also extends into the adjacent provinces which are part of the watershed to ensure that the hydrologic system is fully covered. Administrative boundaries serve as fundamental units to implement water management and availability assessment of water at sub-district level.

2.2. Hydrologic Modeling

2.2.1. QSWAT Model

QSWAT, which integrates Soil and Water Assessment Tool (SWAT) with QGIS, is a comprehensive hydrologic model system designed to simulate hydrologic processes at watershed level [19]. This integrated model was developed to increase capacity in accessing hydrologic analysis, and it has been extensively applied in different fields as well as fluid dynamics assessment, sediment yield analysis, and analysis of impacts of land use and climate change on water resource [20,21,22,28,29,30]. The model depends on hydrologic foundation using water balance equations calculating change in storage as component difference of flow in and out, and it also uses energy balance approach for evapotranspiration modeling. The latest application demonstrates diverse models in different geographic and hydrologic contexts as well as runoff and erosion assessment [44,45], flood risk model development [46,47], efficiency improvement of watershed management [48] emphasizing potential of the model in its accuracy and adaptability. Water balance in QSWAT is shown in Equation (1).

$${\mathit{SW}}_{\mathit{total}}{=\mathit{SW}}_{o }+{{\displaystyle \sum }}_{i=1}^t{(R}_d \text{-} Q_{s }\text{-} E_{a }{\text{-} w}_{\mathit{seep}} \text{-} Q_{\mathit{qw}})$$

(1)

where SW_total is soil water content at time t (mm); SW₀ is initial soil water content (mm); t is time in day unit; R_d is daily cumulative rainfall on day i (mm); Q_s is cumulative runoff on day i (mm); E_a is actual evapotranspiration on day i (mm); w_seep is infiltration on day i (mm); Q_gw is groundwater content on day i (mm).

2.2.2. Data Collection

Necessary inputs for QSWAT obtained from different government organizations include physical features of watershed, meteorological parameters, hydrologic measurement, and operational data from Ubolrat Reservoir (Figure 1) as presented in Table 1.

2.2.3. Sub-Basin Delineation

The study area is located in Maha Sarakham Province covering one part of all Chi and Mun Basins. The Chi Basin has four sub-watersheds including Chi Watershed 3, Chi Watershed 4/1, lower Phong Watershed 2, and Lower Lam Pao Watershed. For sub-watershed division in QSWAT of the Chi watershed, source of the Chi River and Ubolrat Reservoir’s catchment area have to be considered. The Mun Basin also has four sub-watersheds including Lam Siao Yai Watershed 1, Lam Phang Chu Watershed, Lam Plapphla Watershed, and Lam Tao Watershed. All of them are tributaries of the Mun River. In QSWAT, however, division of additional sub-watersheds for the Mun Watershed is unnecessary. Sub-watershed in QSWAT integrates DEM data with available stream network data. This process results in attaining 90 watersheds in Maha Sarakham Province, which is consistent with sub-district administrative boundary as shown in Figure 1.

2.2.4. Importing Rainfall and Meteorological Data

QSWAT is calibrated and validated using rainfall and meteorology data covering divided sub-watersheds as mentioned in 2.1.1. data obtained from the Thai Meteorological Department (TMD) and the Hydro - Informatics Institute (HII) include station 111 rainfall monitoring stations, 34 temperature monitoring station (minimum and maximum), 34 relative humidity monitoring stations, 8 wind speed monitoring stations, and 5 solar radiation monitoring stations. The 14-year daily data are analyzed (2010 - 2023). Spatial distributions of meteorological variables for each sub-watershed are identified by Thiessen polygon method as shown in Figure 2.

2.2.5. Importing Hydrologic Response Units

Generally, all watershed areas are different in land use, soil types, slope, and land management. A watershed area probably has more than one type of land use, or it uses different ways of land management for each plant depending on soil type and land slope [20]. Each part of the watershed has its own hydrologic parameters, such as curve number (CN) based on land use and soil type, which can lead to different hydrologic outcomes [49].

For this study, QSWAT employed a spatial analysis tool in QGIS importing two layers of GIS data, i.e., land use data in 2019 and soil series data in 2015 recorded by the Land Development Department (LDD). These imported data were overlaid in order to gain new data layers. The new layers in each boundary were from the intersection of two-layer data, referred to as Hydrologic Response Units (HRUs). Each HRU can be imported different hydrologic parameters. This calculation can give more details and reliability to the model. This study land use data were imported by land use type based on land use code in QSWAT database, totaling 19 types while the imported soil series types were defined by code in QSWAT database, totaling 41 types. Then the two-layer data spatially overlaid by QSWAT resulted in the attainment of 414 HRUs. The details of land use map and soil series inputs for the model are presented in Figure 3 (a) and (b).

2.3. WEAP Model

Water Evaluation and Planning System (WEAP) model was developed by Stockholm Environment Institute (SEI) in 1988. It serves as an integrating tool for water resource planning and assessment [32]. The model requires basic data regarding hydrology, meteorology, water resources, water demand, and water infrastructure in the study area. WEAP has been globally used because of its capacity in simulating complex system of water resource, analyzing water management policy, and simulating to assess water management options. Model’s flexibility allowed parameters and conditions to be adaptable for assessing the impacts of climate change, land use change, and water infrastructure development in order to analyze priority-based water allocation as well as water scarcity risk assessment by sector [32,33,34,35,36,37,38,39]. WEAP can be integrated with other hydrologic models, such as SWAT, using the result of streamflow simulation as water resource data for water balance analysis and scarcity assessment [33,50,51]. Moreover, compatibility with groundwater model and water quality can lead to comprehensively quantitative and qualitative analysis as well as higher precision and reliability of the results.

WEAP uses four basic equations for water balance modeling and water scarcity analysis [37], including Equation (2): mass balance for water volume change in a system, Equation (3): water balance at the node for storage volume change, Equation (4): water demand calculation inclusive of efficiency and area, and Equation (5): water scarcity specified by difference between demand and supply as follows:

$${\displaystyle \frac{\mathit{dV}}{\mathit{dt}}} = Q_{\mathit{in}} \text{-} Q_{\mathit{out}}$$

(2)

$$S_{(t+1)}{=S}_{(t)}{+I}_{(t) }{+O}_{(t)}$$

(3)

$$W_D{=\mathit{AWR} \text{×} (1\text{-}e}_f{) \text{×} A}_L$$

(4)

$$W_S{=W}_D {\text{-} W}_A$$

(5)

where dV/dt is rate of water volume change rate per time system; Q_in is water volume flowing into system; Q_out is water volume flowing out system; S_(t+1) and S_(t) is storage volume at time t+1 and t respectively; I_(t) is water content flowing in at time t; O_(t) is water content flowing out at time t; W_D is water demand; AWR is rate of water demand per unit of area; e_f is water use efficiency; A_L is water consuming area; W_S is water scarcity

WEAP requires three main types of input: 1) water supply, 2) spatial data, and 3) water demand data (Table 2), including seven components: 1) demand site node, 2) flow requirement node, 3) reservoir node, 4) river, 5) transmission link, 6) return flow, and 7) streamflow gauge. Imported GIS data as the background of model include watershed boundary, river network, water source, streamflow gauge, irrigation area, road, and sub-district boundary as seen in Figure 4 and Figure 5, respectively.

2.4. Model Performance Assessment

Efficiency of QSWAT and WEAP models were assessed by three statistical indices. For the assessment, QSWAT used Coefficient of Determination (R²) to identify the relationship between measured and modeled streamflow values (Equation 6), Nash-Sutcliffe Efficiency (NSE) to measure model’s forecast precision compared to observed mean (Equation 7), and ratio of RMSE to standard deviation of observed values (RSR) (Equation 8). The assessment covered the whole phase of calibration (2010-2016) and validation (2017-2023) - each of them spanned 7 years.

On the other hand, WEAP utilized simulation results obtained from QSWAT as its inputs. The efficiency assessment emphasized R² and NSE values since they have been normally used for precision assessment of streamflow simulation in the water balance model. The results obtained from simulation were compared to the measured data during a 14-year period (2010-2023). The study used statistically reliable criteria for calibrating monthly model to assess its precision and accuracy [52]:

$$R^{2 }= {\displaystyle \frac{{\mathrm{\Sigma }}_{i=1}^n{(O}_i{\text{-}\overline{O})(P}_i\text{-}\overline{P})}{\sqrt{{\mathrm{\Sigma }}_{i=1}^n{{(O}_i\text{-}\overline{O})}^2}\sqrt{{\mathrm{\Sigma }}_{i=1}^n{{(P}_i\text{-}\overline{P})}^2}}}$$

(6)

$$\mathit{NSE}={\displaystyle \frac{{\Sigma }_{i=1}^n{{(O}_i\text{-}\overline{O})}^2\text{-} {\Sigma }_{i=1}^n{{(P}_i\text{-}\overline{P})}^2}{{\Sigma }_{i=1}^n{{(O}_i\text{-}\overline{O})}^2}}$$

(7)

$$\mathit{RSR}={\displaystyle \frac{{\Sigma }_{i=1}^n{{(O}_i{\text{-}P}_i)}^2}{\sqrt{{\Sigma }_{i=1}^n{{(O}_i\text{-}\overline{O})}^2}}}$$

(8)

where O_i is streamflow content recorded at the streamflow gauge; Ō is average streamflow content recorded at the streamflow gauge; P_i is streamflow content calculated by QSWAT and WEAP; $\overline{P}$ is average streamflow content calculated by QSWAT and WEAP; i is data sequence

2.5. QSWAT-WEAP Framework for Water Balance Simulation

According to the concept of water balance assessment by developing integrated QSWAT-WEAP model (see Figure 6-Framework for water balance assessment), the study included both supply and demand sides. The supply side consists of streamflow content simulation of QSWAT that required important inputs, such as DEM, land use, river network, soil type, sub-watershed boundary, meteorological data, reservoir specifications, data regarding storage capacity, storage curve, and evaporation rate, and operational rules. These data were processed in QSWAT to simulate streamflow for each sub-watershed. On the other hand, the demand side, as seen in the upper part of Figure 6, comprises water demand in four sectors, i.e., industry, consumption, agriculture (both irrigated and non-irrigated areas), and services. These sectors were calculated as water demand combined with spatial data including water source, administrative boundary, watershed boundary, irrigation zone, and streamflow gauge location for water balance simulation in WEAP.

The simulation results from WEAP demonstrating on the right side created 4 major spatial outcomes at sub-district level including 1) sub-district streamflow, 2) water demand, 3) water shortage, and 4) level of water shortage. These outcomes were presented through spatial distribution to support decision making on water management and water allocation planning at sub-district level.

3. Results and Discussions

3.1. Streamflow Simulation from QSWAT

3.1.1. Parameter Sensitivity Analysis

Sensitivity analysis of parameters for QSWAT was conducted using SWAT-CUP through Sequential Uncertainty Fitting version 2 (SUFI2) [53] to analyze parameters related to the streamflow calculation in the model compared to 14-year (2010-2023) monthly streamflow data of six monitoring stations. Optimization to obtain comparable results to measured value was 300 iterations by selecting 14 parameters affecting the streamflow [28] and defining Nash-Sutcliffe (NSE) as an objective function. The result of global sensitivity analysis from SWAT-CUP was expressed through t-Stat and p-Value obtained from hypothesis testing of a relationship between parameters and objective function. This analysis indicated that the leading three parameters affecting calculation result were GWQMN (3537), SOL_AWC (0.07), and GW_REVAP (0.83). In analyzing absolute values of t-Stat and p-Value approaching zero, it indicated significant sensitivity to objective function [54]. The remaining 11 parameters with low sensitivity were ranked and presented in Table 3.

3.1.2. Performance Assessment of Models

The result of performance assessment of QSWAT and WEAP models are illustrated in Figure 7, which presents heat map indices, and Figure 8, which represents the consistency of monthly streamflow data from monitoring stations and the calculation results gained from time series model of six stations during 2010-2023. According to the result from QSWAT during calibration period (2010-2016), most stations had R² between high and very high level (0.66–0.87). Station E.66 had the highest level (0.87). NSE was between 0.14 and 0.83; station E22B and Ubolrat Reservoir (Inflow UB) had variation below satisfactory range (at 0.14 and 0.46, respectively). The period of low streamflow was noticed at both stations. The values obtained from the models were higher than the observed values as shown in Figure 8 (b) and 8 (c); however, based on the average of all six stations, NSE was at satisfactory level (0.58). RSR index was consistent with NSE (0.63), which was satisfactory. In the period of model validation (2017-2023), each index showed better outcomes. Stations E9, E91, and E66 indicated their results at very good level from three indices. The remaining stations were at good and satisfactory levels represented by the average of R² = 0.72, NSE = 0.70, and RSR = 0.54.

Performance assessment of streamflow simulation using WEAP was to reverify the calculation results of streamflow gained from QSWAT. In WEAP context, the streamflow of main river, which flows through the study area and become potential and applicable, has led to water balance analysis for three parts: 1) water availability, 2) water demand, and 3) sub-district level water shortage in Maha Sarakham Province. In assessment, R² and NSE indices were used focusing on four main stations (E.9, E.91, E.66A, and M.95). The first three stations aimed to monitor streamflow in the Chi River affecting northern Maha Sarakham Province. On the other hands, M.95 needed performance assessment for Lam Siao River located in the South outside Maha Sarakham Province area. This study found that WEAP was able to simulate streamflow and had R² between 0.66 and 0.80 and NSE between 0.65 and 0.71, compared to four monitoring stations, and had consistency of temporal data during 2010-2023 as presented in Figure 8 (a), (d), (e), and (f). The average of both indices was at good level, and the optimal models ranked by performance were E.9, E.91, E.66, and M.95. The overall performance assessment for QSWAT and WEAP supported that streamflow stimulated by both models provided reliable data from satisfactory to good level compared to model standard and temporal data consistency.

3.2. Spatial Distribution of Water Resources

3.2.1. Sub-District Streamflow Analysis

The results of streamflow analysis in QSWAT and WEAP during 2010-2023 shown in Figure 9 (a) were classified into five levels of spatial distribution. Sub-districts located near the Chi River, which is a main river flowing through northern and central parts of the province (represented by the dark blue in the map), demonstrated high to very high streamflow (Ranges 4-5: 250-7,600 MCM per year) due to surface runoff over sub-district area and influence of water from the main river. On the other hand, southern and southwestern sub-districts in the province (represented by the light green), which is far from the main river and relies only on surface runoff in the area, had low streamflow (Range 1: 2-40 MCM per year). From analyzing quantitative data (see Figure 9 (b)), Range 1 covers 2,212 sq.km or 46% of the total area including 66 sub-districts while Range 5 covers only 359 sq.km or 8% of the total area including 7 sub-districts. For Range 2-4, it covers 2,521 sq.km or 47% including 59 sub-districts (40, 16, and 3 sub-districts respectively) distributing in central part of the province, which is a border area between the main river and areas remote from the river.

The analysis of monthly streamflow changes as shown in Figure 9 (c) and 9 (d) revealed seasonal changes consistent with spatial distribution in Figure 8(a). Ranges 1-3 (Figure 9 (c)) had low streamflow; maximum during September – October with 4.2 - 42.7 MCM streamflow and minimum during February – April with 0.02-0.28 MCM streamflow. This is related to the light yellow, light green, and light blue areas in the distribution map covering 122 sub-districts or 91% of the total sub-districts. In case of Ranges 4-5 as shown in Figure 9 (d), the maximum streamflow was 1,665.5-1,750.3 MCM in October, and the average streamflow in the dry season (December - July) was 94.3-555.6 MCM, which is in line with light and dark blue areas in the map of distribution along the Chi River.

According to the analysis of spatial streamflow distribution corresponding to area size, the number of district, and monthly temporal changes in Maha Sarakham Province, the study found seasonal imbalance and limited streamflow in the dry season. In contrast, intensity became high in the rainy season, which was significantly consistent with spatial contribution and geographical characteristics of Maha Sarakham Province. In the dry season, the areas with 0 – 100 MCM per year streamflow (Ranges 1-2) covering 2,047 sq.km area (78% of the total area) including 106 sub-districts (out of 133 sub-districts) indicated low monthly streamflow that reflected widespread distribution of sub-district level drought. For the central area (100 – 250 MCM per year), it included 709 sq.km or 13% of the total area. The areas with high streamflow (250 - 7,600 MCM per year), which are mostly located near the Chi River covering 489 sq.km including only 11 sub-districts, obviously showed higher value from June to October. This imbalance reflets that the area between central and southern parts had low level of monthly streamflow change while there were not many northern sub-districts having very high intensity of streamflow in the rainy season. Significant difference of streamflow between areas close to and far from the main river reflects a key role of the Chi River in being a major water resource of the province as well as influence of Ubolrat Reservoir to control the streamflow in the dry season. It is consistent with the findings of He et al. [55] indicating that the areas far from the main water resource in the Mekong watershed normally encounter water shortage, and the findings of Cheng et al. [56] highlighting that the distance from a main river is a factor to identify the availability of local water resources.

3.2.2. Sub-District Water Demand Analysis

Water demand at sub-district level through WEAP was analyzed based on five activities: non-irrigated agriculture, irrigated agriculture, consumption, services, and industry as shown in Figure 10 (a) and 10 (b), and it was also divided into five levels. Most of the sub-districts in the province were at level 1 of water demand (1 – 10 MCM per year) including 32 sub-districts and covering 903 sq.km or 17% of the total area, which distributed outside the province to the north and the east (represented by yellow). There were 98 sub-districts having medium level of water demand (Ranges 2 – 4: 10 – 40 MCM per year) covering 4,211 sq.km or 80% distributing in northern, central, and southern part of the province. The highest level of water demand (Range 5: 40 – 55 MCM per year) was found in only three sub-districts covering 178 sq.km or 3% (represented by dark purple in the map). An important notice from comparing between Figure 10 (a) and Figure 9 (a) is that water demand distribution is obviously different from streamflow distribution. In other words, the sub-districts where water demand is high do not need to be close to a main water resource, so this inconsistency becomes an important cause of water shortage for the areas remote to a main river.

Figure 10 (c) presents the change of monthly average water demand which is in line with local crop calendar. All sub-districts have water demand starting from June, which is the period of in-season rice cultivation. For the sub-districts with level 5 of water demand, they need 9.8 MCM water in August while those with Range 1 of water demand require 1.4 MCM water in the same month. During the dry season (January - May), areas of all levels have significantly lower water demand. The sub-districts with Range 1 have average water demand of 0.9 MCM per month. From analyzing the proportion of water demand by sector (Figure 10 (d)), non-irrigated agriculture shows the highest proportion at 1,713.5 MCM per year (88.7%), followed by irrigated agriculture at 185.8 MCM per year (9.6%), consumption at 18.9 MCM per year (0.98%), industry at 12.0 MCM per year (0.6%), and services at 2.6 MCM per year (0.13%). Figure 10 (e) demonstrates the trend of annual water demand, which has continually increased since 1981 from approximately 1,500 MCM to approximately 2,500-3,000 MCM in 2023, representing an average annual increase of 1.7%, caused by agricultural area expansion and population growth. The water demand proportion for non-irrigated agriculture increasing 88.7% reflects that local agriculturists are reliant mainly on rainfall, which is vulnerable to climate change. Considered in conjunction with Figure 9 (a), southern sub-districts in Maha Sarakham Province, where the streamflow is low, have moderate – high water demand (Ranges 2-4). This indicates the imbalance between demand and supply of water that can lead to water shortage.

3.3. Water Shortage Risk Assessment

3.3.1. Sub-District Water Shortage Analysis

Due to calculating the difference between water demand and available streamflow in each sub-district, sub-district level water shortage can be assessed as shown in Figure 11 (a) and classified into five ranges. There were 28 sub-districts having lowest water shortage (Range 1: 0-2 MCM per year) covering 980 sq.km or 19% of the total area. These sub-districts distribute over the north along the Chi River (represented by light yellow) influenced by high volume of streamflow (Ranges 4-5: 250-7,600 MCM per year). Figure 11 (b) presents 79 sub-districts with moderate water shortage (Ranges 2-3: 2-10 MCM per year) covering 3,172 sq.km or 60% of the total area distributing widely around central part of the province. There were 26 sub-districts facing severe water shortage (Ranges 4-5: 10-35 MCM per year) covering 1,140 sq.km or 22% of the total area distributing in the south and the southwest (represented by dark orange and red). These distributions represent the imbalance mentioned in Section 3.2.1 and Section 3.2.2; in other words, the southern area has low streamflow (Ranges 1-2: 2-100 MCM per year), but its water demand is at moderate to high level (Ranges 2-5: 10-55 MCM per year).

Figure 11 (c) and 11 (d) represent monthly change of water shortage obviously different from streamflow and water demand. The sub-districts with Range 1 as presented in Figure 11 (c) have low constant water shortage throughout the year. The maximum is 0.07 MCM in December while the minimum is 0.0 MCM in April, May, and October. The sub-districts with Ranges 2-5 as shown in Figure 11 (d) have maximum water shortage during November – December; 5.63 MCM in November for Range 5 followed by 2.85 MCM for Range 4. All ranges show the minimum in May (0.01-0.08 MCM), which is in the beginning of the rainy season. An interesting finding is that the critical period of water shortage (November - December) is not the highest water demand period (June - August) as described in Section 3.2.2. On the contrary, the streamflow rapidly decreases in post-harvest period while water is demanded for land preparation for the next cycle. This phenomenon indicates that seasonal transition is the critical period of water shortage in agricultural areas [57], and it suggests the importance of water allocation in order of priority during this critical period [58].

3.3.2. Sub-District Level of Water Shortage Analysis

Analyzing proportion of water shortage severity by percentage of unmet water demand receives more insightful answers than considering only absolute water shortage as described in Section 3.3.1 since the sub-districts with low water demand but high proportion of water shortage are probably affected higher than those with high water demand but low proportion of water shortage. Figure 12 (a) demonstrates five levels of severity classification. The sub-districts with high severity cover 1,664 sq.km or 31.4% including 40 sub-districts and distributing in the south and the east of the province (represented by orange). Secondly, sub-districts with medium severity cover 1,397 sq.km or 26.4% including 32 sub-districts and distributing in central part of the province. Figure 12 (b) shows the sub-districts with very high severity. Although it includes only 33 sub-districts, it covers 1,242 sq.km or 23.5% distributing in the southernmost and southwestern parts (represented by red), and overlapping with the areas with low streamflow (Ranges 1-2 in Section 3.2.1). On the other hand, there are only 28 sub-districts with very low to low severity covering 989 sq.km or 16% and distributing in the north along the Chi River.

Monthly severity changes as shown in Figure 12 (c) and 12 (d) reveal a critical situation for the areas with high and very high severity. The sub-districts with very high severity show 99.24% water shortage in March and the severity is still higher than 95% during January – April. It means that water demand has not been responded in the dry season although in the rainy season, the value is still at 35.23 - 51.80%, which is regarded as very high. The sub-districts with high severity show similar result; maximum = 96.82% in March while minimum = 4.96% in October, while those with medium severity have the maximum at 85.00% in March and the minimum at 0.79% in October. Those with very low severity have low and constant values throughout the year, and the maximum is only at 0.22% in June. The overall result indicates that there are 105 sub-districts or 79% of the total sub-districts (medium to very high) significantly affected by water shortage risk in terms of water use activities especially during January – April at all levels (except very low) with higher 50% water shortage. This result suggests that areas from central to southern part need urgent measures, such as minor water source development, promotion of water use efficiency, and cultivation planning based on water availability.

3.4. Synthesis of Water Balance and Shortage Risk Findings

The finding integration between distribution of water resources (Section 3.2) and water shortage risk assessment (Section 3.3) represents an obvious relationship between water supply and demand at sub-district level in Maha Sarakham Province. The study found the imbalance between available streamflow and water demand that was significantly different between the north and the south of the province. The northern area along the Chi River has high streamflow (Ranges 4-5: 250-7,600 MCM per year) covering 11 sub-districts and 489 sq.km (9%), and has very low to low severity of water shortage covering 28 sub-districts and 989 sq.km (18.7%). Conversely, from central to southern area has low streamflow (Ranges 1-2: 2-100 MCM per year) covering 106 sub-districts and 4,094 sq.km (78%), and has high to very high severity covering 73 sub-districts and 2,906 sq.km (54.9%). This difference was caused by the concentration of main water sources in the north while water demand distributed throughout the province; in other words, most sub-districts (101 sub-districts or 76%) had medium to high water demand (Ranges 2-5: 10-55 MCM per year) as shown in Table 4.

Based on temporal analysis, interesting findings on monthly water shortage were revealed. Section 3.3.1 has described that maximum water shortage was found during November – December; that is to say Range 5 area had maximum at 5.63 MCM in November. Section 3.3.2 has indicated that the areas with very high level had 99.24% water shortage in March and were higher than 95% during January – April, which means water demand was not met in the dry season. After considering sub-district characteristics of each risk level as presented in Table 5, sub-districts with very high risk had common characteristics, i.e., being located more than 25 km from a main river, streamflow intensity below 0.5 MCM per sq.km, and dependence on rainwater as a main water source, and the critical period was long from November to June. On the other hand, the sub-districts with very low risk were located near or within five km of the Chi River, and they had the streamflow intensity over 100 MCM per sq.km and benefited from the existing irrigation system. These characteristics indicate that a distance from the main river and irrigation system accessibility are factors for specifying risk level that can be used as criteria to screen target areas for future development of water sources. Overall, the findings support the necessary to develop minor water sources and local irrigation system in the south of the province as well as promotion of water use efficiency for non-irrigated agriculture, which needed high water demand at 88.7%.

4. Conclusions

Maha Sarakham Province is an important agricultural area in northeastern Thailand that has encountered persistent water shortage over the past decade, especially in a long dry season that affects agricultural sector significantly. The Chi River is the main provincial water source; however, distributing water to remote areas is limited while agricultural areas and communities have continuously expanded. As a result, average water demand increases 1.7% per year; in other words, it increased from approximately 1,500 MCM per year in 1981 to approximately 2,500 - 3,000 MCM per year in 2023. This study aimed to develop water shortage risk assessment at sub-districts level through QSWAT and WEAP integration. QSWAT was used for streamflow simulation at sub-watershed level using DEM, land use, soil types, and meteorological data. WEAP was used for water balance analysis and water shortage risk assessment at sub-districts level by analyzing water demand in five sectors: non-irrigated agriculture, irrigated agriculture, consumption, services, and industry.

The research findings revealed that QSWAT and WEAP models were efficient. The average R² and NSE during validation period for QSWAT were 0.72 and 0.70 respectively and for WEAP were 0.74 and 0.67 respectively. According to analyzing water resource distribution, the streamflow became localized in the north of province along the Chi River. There were 11 sub-districts (9%) in this area having high streamflow (Ranges 4-5: 250-7,600 MCM per year) while 106 sub-districts from the center to the south (78%) had low streamflow (Ranges 1-2: 2-100 MCM per year). The average water demand at provincial level was 1,932.8 MCM per year; the maximum proportion was at 88.7% in non-irrigated agricultural sector. In assessing water shortage risk, there were 105 sub-districts (79%) with medium to very high risk faced water shortage that could significantly affect water use activities. The areas with very high risk had water shortage at 99.24% in March and higher than 95% during January – April. The maximum water shortage could be found during November – December; in other words, the areas with Range 5 had 5.63 MCM in November.

In summary, integration between QSWAT and WEAP can assess water shortage risk at sub-district level efficiently, and the results can be presented in an understandable map for local users. Classification of sub-district characteristics at each level of risk can help specify important risk factors, such as a distance from the main river and irrigation system accessibility, which can be used as screening criteria for target areas to develop future water sources. The findings can be utilized by concerned organizations, such as DWR, RID, and Local Administrative Organizations, to develop a plan to manage water resources and define measures of drought solution in the target area, particularly for 73 sub-districts with high to very high risk covering 3,044.16 sq.km (54.79%) from central to southern part of the province. In case of further study, groundwater should be integrated with a model to enhance completeness of water resource assessment and to develop a model to predict the effects of climate change and land use change on long-term water shortage risk. In addition, monitoring and real-time warning systems integrated with a model can improve performance in proactive water management for drought-prone areas as well as possibility studies on development of water distribution from surplus areas along the Chi River to deficit areas in the south.