Next Article in Journal
Estimation of Soil Salt and Ion Contents Based on Hyperspectral Remote Sensing Data: A Case Study of Baidunzi Basin, China
Previous Article in Journal
A Forecast-Skill-Based Dynamic Pre-Storm Level Control for Reservoir Flood-Control Operation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Geological Sciences, Faculty of Science, Chiang Mai University, 239 Huaykaew Rd., Chiang Mai 50200, Thailand
2
Department of Mathematics, Faculty of Science, Chiang Mai University, 239 Huaykaew Rd., Chiang Mai 50200, Thailand
3
Centre of Excellence in Mathematics, CHE, Si Ayutthaya Rd., Bangkok 10400, Thailand
4
Advanced Research Center for Computational Simulation (ARCCoS), Chiang Mai University, 239 Huaykaew Rd., Chiang Mai 50200, Thailand
5
Environmental Science Research Center (ESRC), Chiang Mai University, 239 Huaykaew Rd., Chiang Mai 50200, Thailand
*
Author to whom correspondence should be addressed.
Water 2021, 13(4), 558; https://doi.org/10.3390/w13040558
Submission received: 19 January 2021 / Revised: 15 February 2021 / Accepted: 16 February 2021 / Published: 22 February 2021

Abstract

:
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.

1. 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 Mm3/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, MODFLOW, 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.

2. General Information of the Study Area and Previous Studies

2.1. Hydrogeologic Settings

The unconsolidated aquifers of the Chiang Mai groundwater basin cover an area approximately 2800 km2 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].

2.2. 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.

2.3. 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 steady-state, 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 Mm3/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 Mm3.

3. Materials and Methods

3.1. 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).

3.2. 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].

3.3. 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 Mm3/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.

3.4. 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.
K h + q = 0
K h + q = S s h t
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 km2 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 7 and Figure 8 show a top view of the finite difference grids and 4 model layers (with aquifer materials of the model). Temporal discretization consisted of 1 steady-state (SS) stress period (1 day) and 119 one-month transient (TR) stress periods. Each month of the transient stress periods was divided into five equal timesteps. The simulation started from 31 July 2006 to 1 July 2016. 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).

3.5. Model Calibration

PEST [7,8] was used to assist highly parameterized model calibration based on a pilot-point 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 ( Φ )
Φ = i = 1 N ω i x i x i 2 = i = 1 N ω i r i 2
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 forward-difference perturbation technique as follows.
Φ p j Φ ( p j + Δ p j ) Φ ( p j ) Δ p j
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).
RMSE = 1 N i = 1 N ( h i h i ) 2   and   NRMSE   ( % ) = 1 N i = 1 N h i h i h i 2 × 100
Anderson et al. [22] suggested that the calibrated model is acceptable if RMSE is less than 10% of the targeted heads.

4. Results and Discussion

4.1. 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 Mm3/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 Mm3/y or approximately 60,275 m3/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 illustrates the result of the transient model calibration; the time-variable hydraulic heads from August 2006 to July 2016 of some selected wells. The model could capture the transient behavior of this groundwater system where groundwater level fluctuated seasonally (e.g., G135 and MW988). The final parameter values after calibration are shown in 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, Ss) are illustrated in Figure 12 and Figure 13 indicating that the Chiang Mai aquifer systems are highly heterogeneous. Figure 14 and Figure 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.

4.2. 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.

4.3. 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 Sy) 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 (Sy) varied from 10−2–10−1 whereas specific storage (Ss) 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.

4.4. 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 Mm3 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 Mm3/y. That is, a cumulative pumping rate of 51.2 Mm3/y is the value of safe yield for an average of a 2-m drawdown.

5. 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 m2 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 Mm3/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.

Author Contributions

Conceptualization, S.S.; Data curation, S.T.; Formal analysis, S.T.; Funding acquisition, S.S.; Investigation, S.T.; Resources, S.S.; Software, S.T., M.K., and S.S.; Supervision, S.S.; Validation, S.T., M.K., and S.S.; Visualization, S.T.; Writing—original draft, S.T., and S.S.; Writing—review & editing, S.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research work was partially supported by Chiang Mai University; Thailand Science Research and Innovation (TSRI, formerly the Thailand Research Fund, TRF) (Grant number RDG6230005; MRG4980079); the Centre of Excellence in Mathematics, CHE, Thailand; and Science Achievement Scholarship of Thailand (SAST).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Groundwater monitoring data, as well as simulation input files, are available upon request.

Acknowledgments

This research work was partially supported by Chiang Mai University and the Centre of Excellence in Mathematics (CHE). The authors would also like to acknowledge the Department of Groundwater Resources, Ministry of Natural Resources and Environment, Thailand for providing the long-term groundwater monitoring data for this study.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Department of Groundwater Resources (DGR). Thailand’s Conjunctive Groundwater and Surface Water Use Planning: Northern Thailand; Final Report; Department of Groundwater Resources: Bangkok, Thailand, 2010. [Google Scholar]
  2. Department of Groundwater Resources (DGR). The State of Groundwater Annual Report: Thailand Groundwater Monitoring Networks Project; Final Report; Department of Groundwater Resources: Bangkok, Thailand, 2013. [Google Scholar]
  3. Department of Groundwater Resources (DGR). The State of Groundwater Annual Report: Thailand Groundwater Monitoring Networks Project; Final Report; Department of Groundwater Resources: Bangkok, Thailand, 2018. [Google Scholar]
  4. Department of Groundwater Resources (DGR). Groundwater and Environmental Development and Conservation Master Plan; Final Report; Department of Groundwater Resources: Bangkok, Thailand, 2012. [Google Scholar]
  5. Intasutra, T. Groundwater potential in unconsolidated sediments in Chiang Mai basin. In Proceedings of the First Symposium on Geomorphology and Quaternary Geology of Thailand, Bangkok, Thailand, 28–29 October 1983; pp. 179–195. [Google Scholar]
  6. Department of Groundwater Resources (DGR). Groundwater Resources Assessment of the Chiang Mai, Upper Chao Phraya, and Mae Klong Basins; Final Report; Department of Groundwater Resources: Bangkok, Thailand, 2005. [Google Scholar]
  7. Doherty, J. PEST, Model-Independent Parameter Estimation-User Manual, 5th ed.; Watermark Numerical Computing: Brisbane, Australia, 2014. [Google Scholar]
  8. Doherty, J.; Hunt, R.J. Approaches to Highly Parameterized Inversion: A Guide to Using PEST for Groundwater-Model Calibration: U.S. Geological Survey Scientific Investigations Report 2010-5169; U.S. Geological Survey: Reston, VA, USA, 2010. [Google Scholar]
  9. Margane, A.; Tatong, T. Aspects of the hydrogeology of the Chiang Mai-Lamphun basin, Thailand that are Important for the Groundwater Management. Z. Angew. Geol. 1999, 45, 188–197. [Google Scholar]
  10. Wattananikorn, K.; Beshir, J.A.; Nochaiwong, A. Gravity interpretation of Chiang Mai basin, northern Thailand: Concentrating on Ban Thung Sieo area. J. SE Asian Earth Sci. 1995, 12, 53–64. [Google Scholar] [CrossRef]
  11. Asnachinda, P. Hydrogeochemistry of the Chiang Mai basin, northern Thailand. J. Asian Earth Sci. 1997, 15, 317–326. [Google Scholar] [CrossRef]
  12. Ridd, M.R.; Barber, A.J.; Crow, M.J. The Geology of Thailand; Geological Society: London, UK, 2011. [Google Scholar]
  13. Suvagondha, F. Hydrogeology of Amphoe Muang Lamphun. Master’s Thesis, Chiang Mai University, Chiang Mai, Thailand, 1979. [Google Scholar]
  14. Kiwduangta, B. Hydrostratigraphic Model of the Shallow Aquifer in the Vicinity of the Northern Region Industrial Estate Area, Lamphun Province. Master’s Thesis, Chiang Mai University, Chiang Mai, Thailand, 2012. [Google Scholar]
  15. Wongpornchai, P. Hydrogeology of the Western Part of Chiang Mai Basin: Amphoe Mae Rim, Amphoe Muang and Amphoe Hang Dong, Chiang Mai. Master’s Thesis, Chiang Mai University, Chiang Mai, Thailand, 1990. [Google Scholar]
  16. Tatong, T. Assessment of Available Groundwater Resources in the Chiang Mai Basin, Northern Thailand. Master’s Thesis, University of Birmingham, Birmingham, UK, 2000. [Google Scholar]
  17. Uppasit, S. Groundwater Recharge Calculation of Chiang Mai Basin Using Water-Table Fluctuation Method. Master’s Thesis, Chiang Mai University, Chiang Mai, Thailand, 2004. [Google Scholar]
  18. Saenton, S. Applications of Stochastic Process and Inverse Modeling Technique in Groundwater Modeling of the Chiang Mai Basin, Thailand Research Fund; Final Report; Thailand Research Fund: Bangkok, Thailand, 2010. [Google Scholar]
  19. Harbaugh, A.W. MODFLOW-2005, the U.S. Geological Survey Modular Groundwater Model: The Ground-Water Flow Process; U.S. Geological Survey Techniques and Methods 6-A16; US Department of the Interior: Reston, VA, USA, 2005. [Google Scholar]
  20. Hill, M.C.; Tiedeman, C.R. Effective Groundwater Model Calibration, with Analysis of Sensitivity, Predictions and Uncertainty; Wiley: New York, NY, USA, 2007. [Google Scholar]
  21. Doherty, J. Ground water model calibration using pilot points and regularization. Ground Water 2003, 41, 170–177. [Google Scholar] [CrossRef] [PubMed]
  22. Anderson, M.P.; Woessner, W.W.; Hunt, R.J. Applied Groundwater Modeling: Simulation of Flow and Advective Transport, 2nd ed.; Elsevier Inc.: London, UK, 2015. [Google Scholar]
Figure 1. Digital elevation map and boundary of unconsolidated aquifers Chiang Mai basin (~2800 km2).
Figure 1. Digital elevation map and boundary of unconsolidated aquifers Chiang Mai basin (~2800 km2).
Water 13 00558 g001
Figure 2. Hydrogeologic units of the Chiang Mai basin (modified from Reference [12]), overlain by model grids used in groundwater flow modeling.
Figure 2. Hydrogeologic units of the Chiang Mai basin (modified from Reference [12]), overlain by model grids used in groundwater flow modeling.
Water 13 00558 g002
Figure 3. Hydrogeologic cross-sections of the unconsolidated aquifers (Q’s) in lines AA′, BB′, and CC′. Note that, lower units are grouped into single semi- to consolidated aquifers.
Figure 3. Hydrogeologic cross-sections of the unconsolidated aquifers (Q’s) in lines AA′, BB′, and CC′. Note that, lower units are grouped into single semi- to consolidated aquifers.
Water 13 00558 g003
Figure 4. The conceptual model of the Chiang Mai basin and general groundwater flow directions in Qfd, Qyt + Qot, and lower aquifers (a) and hydraulic head distribution and general groundwater flow directions in aquifer units in the year 2015 (b) and hydraulic conductivity field interpolated from pumping test results and (c) Hydraulic conductivity field from pumping test results.
Figure 4. The conceptual model of the Chiang Mai basin and general groundwater flow directions in Qfd, Qyt + Qot, and lower aquifers (a) and hydraulic head distribution and general groundwater flow directions in aquifer units in the year 2015 (b) and hydraulic conductivity field interpolated from pumping test results and (c) Hydraulic conductivity field from pumping test results.
Water 13 00558 g004
Figure 5. Locations of groundwater monitoring wells and gaging stations of basin’s main rivers.
Figure 5. Locations of groundwater monitoring wells and gaging stations of basin’s main rivers.
Water 13 00558 g005
Figure 6. Locations of groundwater extraction wells in Qfd, Qyt, and Qot + Qcl units and their extraction rates.
Figure 6. Locations of groundwater extraction wells in Qfd, Qyt, and Qot + Qcl units and their extraction rates.
Water 13 00558 g006
Figure 7. Finite-difference grid used to simulate groundwater flow (top view).
Figure 7. Finite-difference grid used to simulate groundwater flow (top view).
Water 13 00558 g007
Figure 8. Distribution of aquifer units in four model layers.
Figure 8. Distribution of aquifer units in four model layers.
Water 13 00558 g008
Figure 9. (a) Grid blocks for rivers and general-head package and (b) locations of regular pilot points used in parameter estimation for spatially variable parameters such as hydraulic conductivities, specific storage, specific yield, recharge, and maximum evapotranspiration rates.
Figure 9. (a) Grid blocks for rivers and general-head package and (b) locations of regular pilot points used in parameter estimation for spatially variable parameters such as hydraulic conductivities, specific storage, specific yield, recharge, and maximum evapotranspiration rates.
Water 13 00558 g009
Figure 10. Steady-state model calibration results: (a) measured vs. model simulated hydraulic heads and (b) hydraulic head (m a.m.s.l.) contours after steady-state model calibration.
Figure 10. Steady-state model calibration results: (a) measured vs. model simulated hydraulic heads and (b) hydraulic head (m a.m.s.l.) contours after steady-state model calibration.
Water 13 00558 g010
Figure 11. Examples of transient hydraulic heads between modeled and observed heads.
Figure 11. Examples of transient hydraulic heads between modeled and observed heads.
Water 13 00558 g011
Figure 12. Spatial distribution of hydraulic conductivities after model calibration.
Figure 12. Spatial distribution of hydraulic conductivities after model calibration.
Water 13 00558 g012
Figure 13. Spatial distribution of specific storage after model calibration.
Figure 13. Spatial distribution of specific storage after model calibration.
Water 13 00558 g013
Figure 14. Spatial distribution of recharge after model calibration.
Figure 14. Spatial distribution of recharge after model calibration.
Water 13 00558 g014
Figure 15. Spatial distribution of evapotranspiration after model calibration.
Figure 15. Spatial distribution of evapotranspiration after model calibration.
Water 13 00558 g015
Figure 16. Variability of recharge over a period of 2007–2015 (from calibration), with a Pearson’s correlation coefficient between precipitation and recharge of 0.79.
Figure 16. Variability of recharge over a period of 2007–2015 (from calibration), with a Pearson’s correlation coefficient between precipitation and recharge of 0.79.
Water 13 00558 g016
Figure 17. Examples of model validation results showing observed (red circles) vs. simulated heads (black dots).
Figure 17. Examples of model validation results showing observed (red circles) vs. simulated heads (black dots).
Water 13 00558 g017
Figure 18. Groundwater recharge and groundwater storage of the unconsolidated aquifers (Qfd, Qyt, Qot) from 2007–2019.
Figure 18. Groundwater recharge and groundwater storage of the unconsolidated aquifers (Qfd, Qyt, Qot) from 2007–2019.
Water 13 00558 g018
Table 1. Hydraulic conductivities and storage coefficients of grouped hydrostratigraphic units from pumping test analyses.
Table 1. Hydraulic conductivities and storage coefficients of grouped hydrostratigraphic units from pumping test analyses.
Hydrostratigraphic UnitsK (m/d)S (-)
Floodplain, Qfd0.015–180.0166–0.299
Low Terrace, Qyt0.035–300.0153–0.204
High Terrace, Qot (+Qcl)0.043–630.0170–0.233
Semi- to consolidated aquifers (below Q’s)0.12–13.51.5 × 10−5–0 × 10−4
Table 2. Initial model parameter values.
Table 2. Initial model parameter values.
Parameter TypesParameter NameInitial Values Note
Hydraulic conductivity b (m/d)QfdHK10.015–18 Water 13 00558 i001Pilot points a
(Layers 1–4)
QytHK20.035–30 m/d
QotHK30.043–63 m/d
Lower aquifersHK40.12–13.5 m/d
Storage coefficient bQfdSS1 (m−1) SY1 (-)5.4 × 10−5–1.2 × 10−4 0.0166–0.2496
QytSS2 (m−1) SY2 (-)3.1 × 10−5–6.7 × 10−4 0.01532–0.1927
QotSS3 (m−1) SY3 (-)1.8 × 10−5–4.2 × 10−4 0.01703–0.2328
Lower aquifersSS4 (m−1)1.5 × 10−5–2.0 × 10−4
Vertical anisotropyAll unitsVANI0.1 -
Recharge rate (mm/y)-RCH28–112 Water 13 00558 i002Pilot points a
(Layer 1)
Max ET rate (mm/y)-EVT88–147
North and south GHB
conductance factors
-GHB0.25 m2/m/d -
River-bed conductance factorsPing riversRIV_10.25 m2/m/d -
Kuang riversRIV_20.25 m2/m/d -
a Pilot points allow parameter values to vary spatially (i.e., heterogeneous or non-uniform); b Fixed pilot points were added to regular pilot points at the locations where pumping tests were conducted.
Table 3. Final parameter values after calibration.
Table 3. Final parameter values after calibration.
Parameter TypesFinal Parameter ValuesParameter Sensitivity
Geometric MeanMin-Max
Hydraulic conductivity (m/d)Qfd0.04690.0051–0.191.07
Qyt0.10500.0187–2.481.78
Qot0.17480.0016–6.771.64
Lower aquifer0.10590.0003–0.420.23
Specific storage (m−1)Qfd0.000261.03 × 10−6–1.60 × 10−41.22 × 10−2
Qyt0.00024.36 × 10−6–1.31 × 10−44.76 × 10−2
Qot0.00011.78 × 10−6–1.82 × 10−42.37 × 10−2
Lower aquifer0.00032.11 × 10−6–2.17 × 10−49.60 × 10−4
Specific yield (-)Qfd0.02943.37 × 10−2–0.203.32 × 10−3
Qyt0.00296.75 × 10−2–0.108.77 × 10−3
Qot0.01128.35 × 10−2–0.114.32 × 10−3
Lower aquifern/a-
Recharge rate (mm/y)-56.50.041–289.620.71
Maximum ET Rate (mm/y)-37.32067.02–188.230.66
Riverbed conductance factor (m2/m/d)Ping0.45170.37–0.510.12
Kuang0.45280.44–0.470.17
General-head conductance factor (m2/m/d)North0.44174.11 × 10−3
South0.33528.32 × 10−3
Table 4. Safe yield evaluation based on targeted drawdown (from current water level).
Table 4. Safe yield evaluation based on targeted drawdown (from current water level).
ScenarioTargeted Average Drawdown
(m)
Safe Yield or Pumping Rate of Qfd, Qyt, and Qot + Qcl Aquifers
(Mm3/y)
% Increase from Current Usage
Base Case-24.0 (2015)-
Case 11.041.673
Case 22.051.2114
Case 33.054.1125
Case 44.054.9129
Case 55.055.4131
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Taweelarp, S.; Khebchareon, M.; Saenton, S. Evaluation of Groundwater Potential and Safe Yield of Heterogeneous Unconsolidated Aquifers in Chiang Mai Basin, Northern Thailand. Water 2021, 13, 558. https://doi.org/10.3390/w13040558

AMA Style

Taweelarp S, Khebchareon M, Saenton S. Evaluation of Groundwater Potential and Safe Yield of Heterogeneous Unconsolidated Aquifers in Chiang Mai Basin, Northern Thailand. Water. 2021; 13(4):558. https://doi.org/10.3390/w13040558

Chicago/Turabian Style

Taweelarp, Sutthipong, Morrakot Khebchareon, and Schradh Saenton. 2021. "Evaluation of Groundwater Potential and Safe Yield of Heterogeneous Unconsolidated Aquifers in Chiang Mai Basin, Northern Thailand" Water 13, no. 4: 558. https://doi.org/10.3390/w13040558

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop