Evaluation of Groundwater Potential and Safe Yield of Heterogeneous Unconsolidated Aquifers in Chiang Mai Basin, Northern Thailand

Chiang Mai basin has an escalating population growth resulting in high demand for water consumption. Lack of surface water supply in most parts of the basin gives rise to the increasing use of groundwater which leads to a continuous decline in groundwater level in the past decades. This study is the first long-term groundwater monitoring and modeling study that aims at developing a transient, regional groundwater flow model of heterogeneous unconsolidated aquifers based on the MODFLOW program. Long-term groundwater monitoring data from 49 piezometers were used in model calibration and validation. The pilot points technique was used to account for the spatial variability of hydrogeologic parameters of heterogeneous aquifers. The simulation results and statistics showed that most sensitive and significant model parameters were spatially variable hydraulic conductivities and recharge rates. The Chiang Mai basin’s unconsolidated aquifers do not have high potential. The water table and/or potentiometric surface in the southeast and southwest areas of Chiang Mai city were continuously decreasing with no sign of recovery indicating critical groundwater condition and careful management must be considered. Safe yield calculation, based on a 2-m average drawdown threshold, suggested that unconsolidated aquifers of the Chiang Mai basin can sustain overall abstraction rates up to 51.2 Mm3/y or approximately 214% of the current extraction rates.


Introduction
Groundwater is an important source of water supply for domestic, agricultural, and industrial use in Thailand. An increase in population, urbanization, industrial, and agricultural growth has led to increasing demand for groundwater resources in many parts of the country where surface water supply is limited. Chiang Mai is one of the fastest-growing cities where groundwater has been utilized extensively to meet the demand of the population, agriculture, and tourism economy [1][2][3]. The Chiang Mai basin (Figure 1) is an intermontane basin covering parts of the Chiang Mai and Lamphun provinces, Northern Thailand. Outside the municipal areas, most populations mainly rely upon the availability of groundwater. The Department of Groundwater Resources (DGR) of Thailand [3] has forecasted an increasing demand at a groundwater usage rate of 4.6 Mm 3 /y, causing a decrease in groundwater storage and associated hydraulic heads. There has been a significant decline in groundwater levels especially in the south (Hang Dong and San Pa Tong districts) and north and northeast (San Kamphaeng and Mae Rim districts), as a result of agricultural and population growth [1][2][3]. Land use in the Chiang Mai basin has also changed recently. Urbanized and agricultural areas have increased due to population growth causing a significant change in groundwater extraction. Moreover, the Royal Thai Government had initiated nationwide groundwater exploration and development projects to mitigate drought, resulting in a drastic increase use of groundwater resources [4].
Few studies were conducted to assess the groundwater potential, as well as the impact of groundwater exploitation of the Chiang Mai basin [5,6]. Groundwater potential, in this work, refers to (1) the ability of the aquifer to store and transmit groundwater, which can be assessed using field estimates of hydraulic conductivities and storage coefficients from pumping tests, (2) the amount of natural recharge that replenishes the groundwater aquifer annually, and (3) the total pumping rate that does not cause a significant decrease in hydraulic heads (i.e., safe-yield). The last two quantities for groundwater recharge assessment, natural recharge and safe-yield, will be obtained from groundwater modeling results. It should be noted that the quality of groundwater, which is beyond the scope of this paper, is not included in groundwater potential assessment.
This paper presents pioneering work that integrates large-scale field survey hydrogeological data to assess groundwater potential, including safe yield using the long-term transient groundwater flow model of unconsolidated aquifers of the Chiang Mai basin. Observed hydraulic head from 49 piezometers during 2006-2016 is used for calibration, whereas the 2016-2019 data is reserved for model validation. The numerical model, MOD-FLOW, will be coupled with the PEST [7,8] computer code to assist model calibration. The model will be used to assess yearly recharge variability for evaluating groundwater potential and to delineate the areas where groundwater usage is critical.

Hydrogeologic Settings
The unconsolidated aquifers of the Chiang Mai groundwater basin cover an area approximately 2800 km 2 ranging from 18 • 30 to 19 • N and from 98 • 45 to 99 • 15 E. The basin is surrounded by mountain ranges on both east and west sides. The main rivers of the basin, the Ping river, generally run from north to south. Thirty-year average annual precipitation is approximately 1200 mm (1989-2018) and the rainy season typically ranges from June to October. The potential evapotranspiration rate normally exceeds precipitation, except for the months between July and September which is the main period of natural groundwater recharge [6,9]. North-south and east-west dimensions of the basin are approximately 70 and 45 km, respectively. The narrowest portion of the basin is approximately 20 km.
There are several geologic units exposed in the area ranging from Precambrian to recent ages (see Figure 2). Unconsolidated Quaternary alluvium deposits lie in an unconformed manner on older rock formations that cover most of the Chiang Mai basin [10][11][12]. The central part of the basin is dominated by the deposition of silts, sands, and gravels transported by the main rivers. Clayey strata are present throughout the sequence. Groundwater wells in this area are relatively shallow (general depth of approximately 20-70 m below ground surface). Figure 3 illustrates examples of hydrogeologic cross-sections in lines AA , BB , and CC ( Figure 2). Note that the semi-to consolidated aquifers are grouped into a single, heterogeneous, unit (dark red color). These cross-sections are generated from a solid model from several lithologic logs and will subsequently be used in the groundwater flow model setup. The Central Alluvial channel is an area of the highest groundwater exploitation. Hydraulic conductivities of these aquifers evaluated from pumping tests were in the range of 0.01 to approximately 20 m/d [6]. Groundwater quality was generally good with total dissolved solids (TDS) less than 250 mg/L. Fluoride and iron concentrations were generally low, especially in the upper part of the aquifer, owing to a high oxygen concentration and flow velocities [1][2][3].
Several studies found that it was extremely difficult to delineate the sedimentary units in the Chiang Mai basin because correlation was almost impossible due to high heterogeneity [9,13,14]. The hydrostratigraphic model showed the association of sedimentary facies with varying depth and that are intercalated. The correlation between adjacent wells was extremely difficult due to the lack of keybeds or recognizable sequences [14].

Groundwater Recharge
A few studies on groundwater recharge in the Chiang Mai basin were conducted [5,13,[15][16][17]. Suvagondha [13] conducted a hydrogeological survey in the Mueang Lamphun district (eastern part of the basin) that included groundwater recharge estimation based on hydrograph analysis. Groundwater recharge was approximately 273.2 mm/y. Intasutra [5] estimated the overall recharge rate of the Chiang Mai basin using flow net analysis and the recharge rates were found to be in the range of 17-143 mm/y. Wongpornchai [15] estimated the groundwater recharge pattern in the southwestern part of the Chiang Mai basin using hydrograph analysis and the recharge was approximately 65 mm/y. Tatong [16] constructed a one-layer groundwater flow model of the Chiang Mai basin using the MODFLOW program and calibrated the parameters. Based on his steady-state model simulation results, groundwater recharge in the colluvial areas was 25 mm/y, whereas the recharge in the central alluvial, alluvial fan and high terrace area was 293 mm/y. Uppasit [17] used the water table fluctuation method (WTF) to quantify groundwater recharge of the Chiang Mai basin and the average recharge rate was 126 mm/y, which was approximately 11% of the total annual precipitation.

Previous Groundwater Modeling Studies
There were a few previous groundwater modeling studies of the entire Chiang Mai basin. Tatong [16] constructed a one-layer steady-state groundwater flow model of the Chiang Mai basin using the MODFLOW program. The model was constructed under the assumption that the aquifers in this area were unconfined. Modeling results suggested that if the abstraction rate were increased to 20% of the abstraction rate in 1996, groundwater level would dramatically decline in the colluvial aquifer. DGR [6] constructed a steadystate, three-dimensional groundwater flow model to evaluate the groundwater potential of the Chiang Mai basin. The model predicted that the annual water budget and safe yield of the basin were approximately 255 and 125 Mm 3 /y, respectively. Saenton [18] further updated the deterministic groundwater flow model developed by DGR [6] and evaluated the uncertainties associated with hydraulic conductivity heterogeneities using stochastic methods to simulate steady-state flow. His study found that the annual water budget of the basin under steady-state conditions was 241 ± 12 Mm 3 .

Conceptual Model Setup: Hydrostratigraphic Units and Hydraulic Properties of the Aquifers
Hydrogeological data was compiled from this study's extensive field survey and borehole investigations and lithologic logs. Five semi-consolidated to unconsolidated deposits of the Chiang Mai basin were grouped into three hydrostratigraphic units: Qfd, Qyt, and Qot. The floodplain deposits (Qfd) aquifer has a thickness of 0-50 m. The upper part of Qfd is unconfined with an average thickness of 30 m, whereas the bottom part of an aquifer is semi-confined to confined. The low terrace deposits (Qyt) aquifer has a thickness of 0-150 m that can be divided into two sub-units: an unconfined aquifer having an approximate thickness of 30 m located in the east and west parts of the floodplain deposits, while the deeper aquifer located beneath the floodplain deposits is confined. The high terrace deposits (Qot) aquifer, which includes Qcl, have an average thickness of 120 m. The unconfined part of this aquifer is approximately 30 m thick, whereas the aquifer beneath the low-terrace deposits is confined. Figure 4 illustrates a conceptual model of the Chiang Mai basin where groundwater in all aquifer units generally flows toward the central part of the basin and discharges to the main rivers in the Qfd unit. According to the flow directions, as shown in Figure 4, it can be deduced that groundwater recharge mainly occurs in Qot and Qyt units. The model domain similar to Figure 2 was cropped at the north and south boundaries and replaced by the general-head boundaries (GHB) in MODFLOW [19] to compensate for lateral inflow/outflow of groundwater. The conceptual model is shown in Figure 3. More than 200 pumping tests were conducted and datasets were analyzed for hydraulic conductivities and storage coefficients, as shown in Table 1. The ranges of K and S were approximately 1-3 orders of magnitude indicating the heterogeneous nature of the aquifers (see examples of hydraulic conductivity distribution in Figure 4). The values K and S in Table 1 will be interpolated into the finite-difference grid used in subsequent simulations as initial parameter values. However, at locations of pumping tests, known K and S values will be set as fixed values (i.e., fixed pilot points).

Hydraulic Head Measurements and River Discharges and Stages Data
Transient hydraulic head data used for the initial conditions and model calibration were obtained from field measurements (2006-2019) from 49 monitoring wells of various depths (see Figure 5). Hydraulic heads from 2006-2016 were used for model calibration, whereas the data from 2016-2019 was reserved for model validation to ensure the applicability of the model. Baseflow data, also used as flow observations in model calibration, was estimated using baseflow separation analyses of hydrograph data obtained from gauging stations of main rivers located at P.1, P.5, P.67, P.73, P.75, and P.81. The use of flow observation in model calibration will increase model reliability and accuracy [20].

Groundwater Wells
There are several active groundwater wells in the Chiang Mai basin which have been used for domestic, agricultural, and industrial purposes (approximately 2500 wells). The locations of these wells are displayed in Figure 6. The extraction rates of these wells were recorded and will be input into the model. It should be noted that the number of wells, as well as extraction rates, vary yearly. In 2019, the number of pumping wells in unconsolidated aquifers was reported as 2568 wells with total extraction rates of 29.8 Mm 3 /y. The number of wells that are screened in each of the Q units (or combination of Q units), as well as total extraction rates in 2019, is shown in Figure 6. Some wells are screened only in a single, whereas the others are screened over two formations. ("n/a" indicates none). This means, most groundwater was extracted from the Qyt and Qot units.

Numerical Model Simulation
After a conceptual model was established, it was converted into a numerical model to simulate groundwater flow and calculate hydraulic heads that vary spatially and temporally. A computer code MODFLOW [19] was used in this study under the aid of graphical user interface software namely Groundwater Modeling Software or GMS ® program (Aquaveo, LLC.). MODFLOW packages used in the simulation included Layer-property flow (LPF), recharge (RCH), Evapotranspiration (EVT), River (RIV), Multi-node well (MNW), and General-head boundary (GHB) packages. Steady-state and transient groundwater flow equations that MODFLOW program approximates the hydraulic head, h, are shown in Equations (1) and (2), respectively.
where K is the hydraulic conductivity tensor, q refers to the internal source/sink.
The model has four materials: the flood plain unconsolidated sediments aquifer unit (Qfd), young terrace unconsolidated sediments aquifer unit (Qyt), high terrace unconsolidated and semi-consolidated sediment (Qot), and a semi-to consolidated aquifer unit (i.e., lumped lower aquifers). A uniform grid size of 1000 × 1000 square meters consisted of 71 columns, 80 rows, and 4 layers of variable thickness. A uniform grid size of 1 × 1 km 2 was selected because it satisfied numerical convergence criteria and reasonable simulation time. These criteria would ensure that the model calibration step, which involved several hundred to thousands of model runs, would be finished within a reasonable time and would not be crashed or halted during inversion simulation.  Figure 9a shows RIV and GHB grid cells. Table 2 shows the initial parameters used in the model's initial run before calibration.     The RIV package required the input of river length passing through each grid and this is calculated automatically from an input of a digitized arc (.shp file) by GMS. Measured river stages and widths at gaging stations were input in the model using the RIV package in MODFLOW and they were interpolated into finite difference grids. The riverbed (and GHB) conductance was evaluated from model calibration. A multi-node well (MNW) package was used to simulate groundwater extraction, especially from wells that had multiple screen intervals or wells that were screened more than one aquifer unit (or more than one model layer).

Model Calibration
PEST [7,8] was used to assist highly parameterized model calibration based on a pilotpoint technique to capture aquifer heterogeneity, as well as spatially variable recharge and evapotranspiration. The pilot points technique is a method for estimating spatially variable property (such as aquifer properties, recharge, and evapotranspiration). Pilot points are used to substitute zones of parameter uniformity, allowing the disposition of areas of high and low hydraulic property contrast values to be inferred through the calibration process, without the need to define zones before estimating the parameters [21]. A number of the pilot points were introduced to the model domain (see Figure 9b) and PEST was asked to estimate the parameter values of the aquifer at each point. These point values would then be spatially interpolated to all of the active cells within the model domain using an interpolating algorithm allowing heterogeneity of the model. For locations of known hydraulic parameters (from pumping test), the pilot points are forced to remain unchanged and they are called "fixed pilot points." During model calibration, PEST iteratively updates parameter values to minimize the sum of weighted squared residuals (Φ) where x i , ∼ x i and ω i are measured data (heads or flows), simulated equivalents, and weights of observation i, respectively. r i = ∼ x i − x i is residual for observation i. The magnitude of changes in parameter values during each iteration depends on the sensitivity of the parameter. Quantitatively, the sensitivity of a parameter p j , or ∂Φ/∂p j , which indicates the relative importance of the parameter to model output, can be approximated by forwarddifference perturbation technique as follows.
PEST will report parameters' sensitivities after each iteration, as well as at the end of calibration processes.
A goodness of fit of hydraulic heads from model calibration can be assessed using RMSE (root mean squared error) and NRMSE (normalized root mean squared error) values which can be calculated using Equation (5).
Anderson et al. [22] suggested that the calibrated model is acceptable if RMSE is less than 10% of the targeted heads.

Steady-State and Transient Model Calibration
Steady-state groundwater flow model was simulated under equilibrium conditions, such as representing long-term average hydraulic balance, and/or conditions where aquifer storage changes are not significant [22]. The model was also used to investigate hydraulic properties, boundary conditions, long-term (or average) behavior of evapotranspiration and recharge rates, and sensitivity analysis of different model parameters. Figure 10a shows a scatter diagram for comparison of observed (measured in August 2006) versus simulated hydraulic heads. The error bars shown on the plot indicate the minimum and maximum possible heads observed during 2006-2016. Although some wells showed large discrepancies between modeled or simulated heads vs. measure heads (e.g., G265 and MW187), this regional model appears to be relatively reliable with the majority of observation wells fall within a 95% confidence interval limit. The values RMSE and NRMSE are 0.86 m and 0.42%, respectively. Figure 10b shows a steady-state hydraulic head distribution, indicating that groundwater generally flows from the north, east, and west toward the central plain and southward with the major cone of depression located at southern parts of the basin.
Water budget calculation from steady-state simulation indicated that the major source of groundwater was recharge of approximately 140.69 Mm 3 /y, which corresponded to an average recharge of 134 mm/y over the entire basin. The main groundwater outflows were due to evapotranspiration, discharge to main rivers. Groundwater usage through pumping wells in unconsolidated aquifers was 22 Mm 3 /y or approximately 60,275 m 3 /d. River leakage had outflow more than inflow indicating the majority portion of rivers were gaining streams. General-head boundary inflow (in the north) had a high influx of groundwater, as expected because the northern part of the basin is mainly a groundwater recharge area.
The calibrated steady-state model was subsequently used as a starting point (i.e., initial conditions) of a long-term (10-year) transient simulation. Figure 11 Table 3. Spatially variable parameters such as hydraulic conductivities, recharge, and evapotranspiration rates include, geometric mean, minimum and maximum values, and their values were within reasonable range compared with initial parameters estimated from field tests and water budget calculations from previous studies. The true values of the riverbed and general-head conductance, on the other hand, were not available for comparison. Sensitivity analysis indicated the most sensitive parameters were hydraulic conductivities, recharge rates, evapotranspiration, and riverbed conductance, respectively (see Table 3).   Final parameters distributions (K, S s ) are illustrated in Figures 12 and 13 indicating that the Chiang Mai aquifer systems are highly heterogeneous. Figures 14 and 15 illustrate the recharge and evapotranspiration, which vary temporally from year to year, respectively. Zones of high recharge rates are located in the northern, eastern, and parts of southwestern areas of the basin. These areas should be reserved for recharge protection areas. The transient model was investigated further to evaluate the effect of yearly precipitation variation on groundwater natural recharge. Based on the water budget analyses, it was confirmed that the values of natural recharge vary temporally and have a strong dependency on the precipitation pattern (see Figure 16), with a Pearson's correlation coefficient of 0.79. In recent years, the overall precipitation has been declining significantly and so has the groundwater recharge. This is also confirmed by observed groundwater heads at several locations where groundwater table or potentiometric surface, though seasonally varying, are continuously lowered.

Model Validation
After the transient calibration was completed, the model was executed in a forward mode to demonstrate its validity and applicability in prediction using an independent set of hydraulic head data (2016-2019). In the validation simulation, actual pumping rates and river stages were input. The recharge rates during 2016-2019, on the other hand, were approximated from annual precipitation (based on the recharge trend shown in Figure 16, which is approximately 11-12% of the annual precipitation). As shown in Figure 17, the model was able to capture the trends of measured hydraulic heads, which decreased over time.

Groundwater Potential of the Unconsolidated Aquifers
The groundwater potential of the basin has several definitions depending on the methods of evaluation. Traditionally, the groundwater potential of an aquifer system is assessed based on hydraulic parameters such as hydraulic conductivity (K) and storage coefficient (S) obtained from pumping test analyses. Aquifers with high K and S can be said to have high potential because of high flow and storage potentials. Another aspect of groundwater potential can be referred to as the ability of the aquifer to recover from hydraulic head drops. If the water table or potentiometric surface from the observation well(s) rebounds (or recovers) quickly within a reasonable time (e.g., dry and wet seasons or during pumping scheme), such an aquifer can be regarded as high potential. A more quantitative assessment of groundwater potential refers to an estimation of groundwater recharge and groundwater storage (or 'change' in groundwater storage). Total or net recharge estimation of a basin has been conducted using several methods such as groundwater level fluctuations, chloride balance, flow net analysis, water budget, etc., although the spatial and temporal distribution of recharge is not easy to evaluate. Analysis on groundwater storage or 'change' in groundwater storage, on the other hand, cannot be assessed straightforwardly and, thus, a groundwater model can be applied.
The groundwater model does not only provide the aquifer's hydraulic parameters, but it also gives quantitative information on net recharge and aquifer storage. In this study, groundwater potential can be assessed as follows.
Hydraulic conductivity is one of the most important hydrogeological parameters needed for managing groundwater resources. It describes an ability of an aquifer to transmit water per unit area. Although a high value of hydraulic conductivity alone may not guarantee a large quantity of extractable groundwater, it can be used to assess how fast groundwater can flow and recover after pumping. The average hydraulic conductivities of the unconsolidated aquifers (Qfd, Qyt, Qot) of model layers 1-3 were in the range of 10 −3 to 10 −1 m/d (three orders of magnitude), indicating poor to fair aquifer conditions. Storage coefficient is the volume of water that the aquifer releases or takes into storage per unit surface area of the aquifer per unit component of the head normal to that surface. In an unconfined aquifer, the storage coefficient corresponds to its specific yield. The value of the storage coefficient for typical water table aquifers (i.e., specific yield or S y ) varies from 0.01 to 0.35, while for well-confined aquifers, specific storage varies from 10 −5 to 10 −3 [22]. In the Chiang Mai basin, field studies and model calibration results show that specific yield (S y ) varied from 10 −2 -10 −1 whereas specific storage (S s ) shows a wider range of 10 −6 -10 −4 m −1 . Values indicate that the storage coefficient of the Qfd, Qyt, and Qot aquifers of the Chiang Mai basin is generally lower than the normal range, so the groundwater potential is probably not in good condition.
Groundwater recharge and storage: Figure 18 illustrates transient groundwater budget, recharge, and storage during the simulation period. It is clear that, in recent years, groundwater recharge and storage of the unconsolidated aquifers (Q's) is declining. This is partly due to the lowering of annual precipitation, as discussed earlier. Based on the above assessment, it can be deduced that the groundwater potential of the unconsolidated aquifers (Qfd, Qyt, Qot), where the average total thickness is~125 m, is not in good condition. Careful groundwater usage must be taken into consideration because groundwater storage is declining according to increasing groundwater usage and decreasing recharge. This observation is confirmed by a continuous decline in the hydraulic head which is indicated in long-term groundwater monitoring over the past fifteen years.

Evaluating Groundwater Safe Yield in Unconsolidated Aquifers
Sustainable or safe yield refers to the amount of groundwater that can be extracted from an aquifer on a sustained basis, economically and legally, without deteriorating the native groundwater quality or creating an undesirable effect, such as environmental damage. In this work, the safe yield is defined as the amount of extractable groundwater through pumping wells that may cause the drawdown, on average, to drop not more than the specified values. The pumping rate from Q units of 24.0 Mm 3 in 2015 was used as a base case. Then, the pumping rates were incrementally increased until the drawdown in all active extraction wells dropped, on average, to the specified preset criteria. Table 4 shows that the values of yield (i.e., groundwater abstraction or pumping rates) corresponded to the targeted average drawdown of the basin for 1-, 2-, 3-, 4-and 5-m drawdown. For example, if the targeted drawdown was set at 2-m below current heads, groundwater abstraction should increase from 24.0 to 51.2 Mm 3 /y. That is, a cumulative pumping rate of 51.2 Mm 3 /y is the value of safe yield for an average of a 2-m drawdown.

Conclusions
This work represents a long-term, up-to-date, and comprehensive groundwater modeling study of the Chiang Mai basin. The model was set up in the simplest, most effective, manner because of data limitations, especially on the availability of hydraulic conductivity distribution and actual groundwater wells. A transient finite-difference groundwater flow simulation of the Chiang Mai basin was setup using the MODFLOW program. The model primarily focused on unconsolidated aquifers (Qfd, Qyt, and Qot), where the majority of local or domestic groundwater wells are extracted. The model calibration was achieved based on an automated parameter estimation scheme using the PEST program. Parameter sensitivity analysis during calibration was obtained using the perturbation method. The sustainable or safe yield based on specified average drawdown, as well as the water budget or groundwater potential, was assessed.
The numerical model consisted of 71 columns, 80 rows, and 4 layers. A uniform horizontal grid size of 1000 × 1000 m 2 with variable thickness for all model layers was used. Major aquifers of this basin are unconsolidated to semi-consolidated aquifers. Model simulation results indicated that hydraulic heads in the north, east, and west regions were high and continuously decreased toward the center of the basin, and groundwater was discharged to main rivers and flowed southward. In some areas, groundwater level dropped significantly below ground surface, especially in the central plain where groundwater abstraction was high due to agricultural, domestic, and industrial needs. The model was calibrated using an automatic parameter estimation program namely PEST with the use of the pilot point technique to evaluate heterogeneous or spatially variable hydraulic conductivities, recharge rates, maximum evapotranspiration rates. Parameters' sensitivity analysis from PEST allowed the assessment of the relative importance of the model parameters. The results showed that the model outcome was sensitive to, in order of importance, hydraulic conductivities, recharge rates, and evapotranspiration rates. The root mean square (RMSE) and normalized root-mean-square (NRMSE) errors of the steady-state model, based on 49 observation wells, were 1.26 m and 0.62%, respectively.
The transient model was calibrated and validated to demonstrate the capability of the model to simulate the seasonal variation of the hydraulic head (2006-2019). The model also indicated the prior assumption where recharge was expected to occur from July to September. Based on the water budget analysis from the transient model, it was confirmed that the values of natural recharge vary temporally and has a strong dependency on precipitation pattern. In recent years, the overall precipitation is significantly declining and so has the groundwater recharge. This is also confirmed by observed groundwater heads at several locations, where groundwater table or potentiometric surface, though seasonally varying, are continuously lowered.
Safe yield calculation suggested that the basin can sustain abstraction rates up to 214% (~51.2 Mm 3 /y for Qfd, Qyt, Qot) of the current extraction rates for a 2-m drawdown threshold. Overall, the groundwater potential of unconsolidated aquifers (Qfd, Qyt, Qot) in the Chiang Mai basin is generally not in good condition, indicating a crisis may occur if uncontrolled use of groundwater of the basin is not properly managed.
Although the numerical model of the regional groundwater system presented in this study was satisfactorily calibrated and validated and able to capture the transient head trends, it was clear that the model did not produce a perfect match in some wells where discrepancies between observed and simulated heads still exist. The model should then be used with caution. While the average or overall response of this regional groundwater system such as storage and recharge can be evaluated, there are still uncertainties associated with heads in individual wells, and interpretation or application must be handled with care.  Data Availability Statement: Groundwater monitoring data, as well as simulation input files, are available upon request.