Distributions of Groundwater Age under Climate Change of Thailand’s Lower Chao Phraya Basin

Abstract

Groundwater is important for daily life, because it is the largest freshwater source for domestic use and industrial consumption. Sustainable groundwater depends on many parameters: climate change is one factor, which leads to floods and droughts. Distribution of groundwater age indicates groundwater velocity, recharge rate and risk assessment. We developed transient 3D mathematical models, i.e., MODFLOW and MODPATH, to measure the distributions of groundwater age, impacted by climate change (IPSL-CM5A-MR), based on representative concentration pathways, defined in terms of atmospheric CO₂ concentration, e.g., 2.6 to 8.5, for the periods 2020 to 2099. The distributions of groundwater age varied from 100 to 100,000 years, with the mean groundwater age ~11,000 years, generated by climate led change in recharge to and pumping from the groundwater. Interestingly, under increasing recharge scenarios, the mean age, in the groundwater age distribution, decreased slightly in the shallow aquifers, but increased in deep aquifers, indicating that the new water was in shallow aquifers. On the other hand, under decreasing recharge scenarios, groundwater age increased significantly, both shallow and deep aquifers, because the decrease in recharge caused longer residence times and lower velocity flows. However, the overall mean groundwater age gradually increased, because the groundwater mixed in both shallow and deep aquifers. Decreased recharge, in simulation, led to increased groundwater age; thus groundwater may become a nonrenewable groundwater. Nonrenewable groundwater should be carefully managed, because, if old groundwater is pumped, it cannot be restored, with a detriment to human life.

1. Introduction

Water plays an important role for life, especially groundwater—the largest freshwater source. It is a source of freshwater for human consumption, agriculture and industry, especially during periods of critical drought and disaster [1]. However, climate change impacts every region, such as agriculture and water resources [2,3,4]. In developing countries, groundwater may be facing shortages and contamination, because there is no plan to manage groundwater under climate change [5,6].

Climate change has an impact on a global scale [7,8]. The change in temperature and rainfall pattern are controlled by climate change [9,10]. In addition, climate change has caused fluctuations in the hydrologic cycle, which impacted the sustainability of surface and groundwater. In particular, tropical climate changes strongly affect groundwater quality and quantity [10,11,12,13], due to the strongly alternating hot and rainy conditions, indicating that climate change directly impacts both temperature and rainfall, which, in turn, affects recharge and hydraulic head [14,15]. The change in groundwater resource under climate change has been studied; for example, Jackson et al. [16] used the A₂ emission scenario of the Intergovernmental Panel on Climate Change (IPCC) General Circulation Models (GCMs), which assumed a population growth of 15 billion by 2100 and a significant decline in fertility for most regions, but stabilizing above replacement levels [17], to project climate by simulating the Chalk aquifer in the UK; they applied a distributed recharge in groundwater flow model and concluded that future groundwater recharge would range between −26% and 31% until 2080. Doll [18] determined that groundwater recharge, under climate change in Northern Brazil, the southern edge of the Mediterranean Sea and southwest Africa, was decreased to 70% of groundwater recharge, under the A₂ and B₂ emission scenarios by 2050. Emelyanova et al. [19] showed that the groundwater levels slightly rose in Southwestern Australia. Shrestha et al. [20] described groundwater levels in the Mekong Delta, in Southern Vietnam, under climate change, which showed a declining trend. These studies showed that change in groundwater recharge and groundwater level depended on the selected emission scenarios and events [21].

Climate change is a very important factor in projecting groundwater recharge [22,23]. However, groundwater sensitivity to climate change is influenced with the aquifer type and depth of groundwater, with both shallow and deep, and small and large aquifers, varying with climate change [24]. Then, the recharge and hydraulic head analysis was not sufficient to understand the impact of climate change on groundwater. The new indicator is groundwater age, because it can be analyzed by groundwater travel time and is not associated with aquifer type, depth and size.

Groundwater age, i.e., the time between water infilled to the aquifer until that water reached a position of interest [25,26], has been shown to be a reasonable indicator of groundwater recharge [27] and sensitivity of groundwater to the climate change and has been assessed in several complex systems [28,29,30,31,32,33,34]. Torgersen et al. described the ideal groundwater age as the time elapsed between when water entered the saturated zone to the time the water was sampled at a specific location, i.e., a specific distance downstream, and defined old groundwater as groundwater that was recharged on a timescale from approximately 1000 to more than 1,000,000 years [35]. The first approach to groundwater age modelling was described by Goode [36], who opened many aspects in the field of numerical groundwater research. Initially, he used the mean over one basin, i.e., only one age. Ginn [37] later showed that groundwater age was distributed over a region and transported in many directions, so the groundwater age is the age assigned to each particle of water. Vani and Carrera [38] showed that when purely advection was considered, the age distribution increased drastically, because the age distribution needed a very fine grid to model the transport: they obtained a more realistic aquifer using a very high resolution (i.e., space, time and age dimensions). Woolfenden and Ginn [39] modeled a real-world groundwater age distribution and showed the modeled age, measured to confirm the groundwater calibration, led to confidence in groundwater age distributions [37] and mean ages [36], calculated by the model.

Suckow [40] reviewed groundwater age dating and showed the advantage of using it to estimate recharge rate, flow rate and model calibration. Groundwater age can indicate groundwater recharge, i.e., young groundwater indicates high recharge, so groundwater can be pumped, but old groundwater indicates low recharge; pumping should be carefully managed because groundwater may not flow to the pumping location. Groundwater age can indicate groundwater flow direction, i.e., if we know the groundwater age, in many locations, we may plot groundwater contours. Groundwater flows from young to old age. Groundwater age can indicate sustainability of groundwater, i.e., if groundwater moves from a young and to an old age, groundwater will become inadequate, due to slower recharge. Then, groundwater should be conserved. In this research, groundwater age was applied to indicate the sustainability of groundwater.

Groundwater age was applied, in many numerical groundwater studies, for example, to estimate groundwater parameters. Sanford used groundwater age to improve overall calibration, such as variability in recharge, porosity and hydraulic conductivity [41]. Michael and Voss used groundwater age to estimate the anisotropy, using MODFLOW and MODPATH in inverse modeling, and found that increasing hydraulic conductivity led to decreased age at each point, and found that in the zone of a discontinuous aquitard, or gap in an aquitard, groundwater was young [42]. The groundwater age may also be combined with other parameters for calibrating groundwater modeling. Portniaguine and Solomon used groundwater age and head data to calibrate a groundwater model, because at some locations, where the head data were not available, the groundwater age streamlines were used to imply to groundwater flux, from which the recharge and hydraulic conductivity may be derived [43]. Engdahl et al. used groundwater age distributions to estimate effectively homogeneous and heterogeneous aquifer: in a homogeneous aquifer, estimation of the properties at all scales was generally accurate, but in heterogeneous aquifers, it only provided reasonable orders of magnitude [44].

Suckow [36] reviewed methods to measure groundwater age; for example, radioactive isotope tracking [35,45] and groundwater models [39,46,47,48]. The methods to measure groundwater were applied to aquifer recharge. Engdahl and Maxwell [49] measured the change in age distribution under variable recharge (i.e., ±25% and −50%) and showed that the change in recharge and topography controlled the age distribution in steady-state conditions, because the recharge rate contributed to groundwater flux (velocity): low flux led to increased groundwater age, whereas high flux contributed to decreased age. Then, “old” groundwater (nonrenewable groundwater) indicated groundwater shortage, when groundwater was pumped, due to a low recharge rate, and it could not be recharged in human lifespans [50].

However, in recent research, it was shown that climate change impacted surface water and groundwater, using mathematical model, for example, sensitivity analysis on climate, hydraulic head, water quantity and quality [2,10,12,13,51]. The major parameters were rainfall and temperature. However, no one focused on the climate change impact on the distribution of groundwater age, which is the contribution of this study.

Thus, we developed a three-dimensional groundwater flow model to investigate the groundwater shortage affected by climate change, using groundwater age as an indicator. First, we describe the study location, including topography, hydrogeology and climate. Then, we detailed the materials and method for calibration and verification of the groundwater model. Since groundwater recharge predictions depend on the selected groundwater simulation model, discussion of the recharge rate in variable scenarios was included. Then we determined the groundwater age distribution to guide groundwater management under climate change.

2. Materials and Data

2.1. Modelling Framework

A MODular three-dimensional finite-difference groundwater FLOW model (MODFLOW) is a numerical modelling for groundwater systems [52]. MODFLOW solves the governing groundwater equations for three-dimensional saturated flow, combining Darcy’s Law and the conservation of mass principle [53]. The general groundwater flow equation form over the physical dimension is:

$$\frac{\partial }{\partial x}\left[K_x\frac{\partial h}{\partial x}\right]+\frac{\partial }{\partial y}\left[K_y\frac{\partial h}{\partial y}\right]+\frac{\partial }{\partial z}\left[K_z\frac{\partial h}{\partial z}\right]+W=S_s\frac{\partial h}{\partial t}$$

(1)

where K_x, K_y and K_z are the hydraulic conductivities in x, y, z directions (LT⁻¹), h is the hydraulic head (L), W is sink or source (T⁻¹) and S_s is specific storage (L⁻¹). Equation (1) shows that MODFLOW can simulate external sources; for example, pumping, recharge, evapotranspiration, drains and rivers. However, MODFLOW does not include the distribution of groundwater age. Then, a Particle-Tracking Model (MODPATH) [54] was adapted to include this factor; the solution of the particle tracking equation over the physical dimension is

$${\left(x_p\right)}_{t_2}=x_1+\frac{1}{A_x}\left[{\left(v_{x_p}\right)}_{t_1}{e}^{\left(A_x∆t-v_{x_1}\right)}\right]$$

(2)

$${\left(y_p\right)}_{t_2}=y_1+\frac{1}{A_y}\left[{\left(v_{y_p}\right)}_{t_1}{e}^{\left(A_y∆t-v_{y_1}\right)}\right]$$

(3)

$${\left(z_p\right)}_{t_2}=z_1+\frac{1}{A_z}\left[{\left(v_{z_p}\right)}_{t_1}{e}^{\left(A_z∆t-v_{z_1}\right)}\right]$$

(4)

where x_p, y_p, z_p are the last particle positions (L), x₁, y₁, z₁ are the initial positions (L), A_x, A_y, A_z are velocity gradients, v_x, v_y, v_z are velocity in x, y, z dimensions, ∆t is the time between the last and initial particle positions, which indicate groundwater age.

MODPATH was used for variable groundwater conditions. Fioreze and Mancuso [55] used MODFLOW-MODPATH to check the groundwater flow pattern in horizontal subsurface flow constructed wetlands. Wolfgang-Albert and Michl [56] applied MODFLOW-MODPATH to check the transport of groundwater, under boundary conditions, e.g., river and recharge, in the alluvial aquifers in Germany and showed the good matching of groundwater levels. Mondal and Singh [57] used MODFLOW-MODPATH to predict contaminant transport in groundwater in India and revealed that, although the pollution levels were reduced by 50%, total dissolved solids (TDS) level in the groundwater, even after 20 years, was not reduced to below 50%. Abrams [58] corrected the groundwater model, using MODFLOW-MODPATH, and revealed that accurate modeling of high-resolution models may not be practical for large watersheds or a simple coarse model. Also, Sanford and Buapeng [45] used MODFLOW-MODPATH and ¹⁴C age-dating to assess groundwater transport and showed that the simulated age could be forced to match the observed age by increasing effective porosity.

Although MODPATH is widely used, it still has three limitations [54]:

• The semi-analytical particle-tracking method is available only for a linear velocity interpolation. MODPATH calculated velocity using interpolated velocities from inter-cell flow rates for the finite-difference approximation in the governing equation.
• The path line analysis depended on discretization of a finite-difference.
• The most important limitation is the uncertainty in boundary conditions and hydrogeologic properties. MODPATH analysis uses only information on ideal water movements, derived from MODFLOW.

These limitations were handled by using a finer grid [59], so that the velocity was matched to a higher dimensional curve. The key limitation is the uncertainty in boundary condition and hydrogeologic properties; here we considered the sensitivity of recharge, pumping and constant head boundaries and also considered hydraulic conductivity and anisotropy. See Tanachaichoksirikun et al. [60] to model precise groundwater age behaviors.

Here, we used Modular Finite-Difference Groundwater Flow Model-2000 (MODFLOW-2000) [52] to simulate the groundwater flow field, because the Lower Chao Phraya basin is a complex groundwater system. In addition, Particle-Tracking Model (MODPATH version-6) [54] was used to determine movement and age of groundwater, because it provides information about water movements and is effective for determining changes in groundwater age with MODFLOW simulation.

2.2. Study Area

Thailand’s Lower Chao Phraya (LCP) basin is in the central plain. The area is 43,300 km²: the east–west edge is ~230 km and the north–south edge ~300 km (Figure 1a). The west side of the basin is connected to the Tenasserim hills, a long range of hills; groundwater flows through the Upper Chao Phraya basin in the north side; the east side is bounded by a series of small hills and the Gulf of Thailand in the south side. Additionally, this basin has four large rivers, the Chao Phraya, Pa Sak, Tha Chin and Mae Klong Rivers and the surface water flows from north to south.

The hydrogeology was deposited in the Upper Tertiary to Quaternary of the Lower Central Plains of Thailand (Figure 1b). Sands, gravels and clays of the Pliocene-Pleistocene-Holocene were composed and sedimented into eight principal confined aquifers, separated by confining clay lenses [61,62]. These aquifers were covered by Bangkok Clay. The total depth is ~500 m. The primary targets for water well construction were the 2nd (Phra Pradeang; PD), 3rd (Nakorn Luang; NL), 4th (Nonthaburi; NB) aquifers, with depths ranging from 100 to 250 m below the ground surface [63].

The LCP basin was studied here, because it is the largest groundwater supply for the central part of Thailand and covers Bangkok, the capital city. There are many groundwater monitoring stations, which can be used for calibration and verification. However, this region also has a capital city and surrounding urban areas, so there is considerable industrial, domestic and agricultural water consumption. So, groundwater may be in short supply and contaminated. In addition, this region suffered from excessive groundwater extraction, which contributed to land subsidence [63]. Additionally, the basin is located in a tropical climate region and is thus sensitive to climate change [10,11,12,13]. Thus, the LCP basin is important for studying climate change impacts on groundwater resources: the distribution of groundwater age is an indicator, so we investigated the groundwater response and sensitivity.

2.3. Climate

The climate in the study area is tropical: the average daily temperature increased from 27.1 °C in 1984–1998 to 27.4 °C in 1999–2014. December is the coolest month with a daily temperature of ~16.6 °C and ~38.3 °C in April, which was the hottest month during 1984–2014 (Figure 2). The average annual rainfall was calculated using ten hydrometeorological stations between 1984 and 2014 (Figure 1a). Almost 60% of rainfall occurred in the rainy season from May to October. The annual rainfall varied from 1000 to 1675 mm (average = 1155 mm). Climate change led to changes in temperature and precipitation patterns [9,10] and, in the tropical LCP basin [10,11,12,13], the recharge rate changed with the climate, which, in turn, led to groundwater age change.

However, in this research, temperature was omitted, because groundwater aquifer was at a 50–100 m depth, thus change in air temperature did not impact the temperature of groundwater.

2.4. Groundwater Age Collection

Sanford and Buapeng [45] measured isotopic ratios for ¹⁸O, ²H, ³H, ¹³C and ¹⁴C, in Thailand, from 1987 to 1996. More than 150 groundwater samples were collected with 1 m well screen lengths, so samples were taken from one aquifer only. About 50 samples were analyzed for ¹⁴C, which is the indicator for the groundwater age, because ¹⁴C was not contaminated by either ³H or seawater and the regional groundwater is in the 5730-year range for the ¹⁴C half-life. Groundwater age from ¹⁴C dating was calculated from:

$$Age=\frac{5730}{ln2}ln\left(\frac{100A_0}{A}\right)$$

(5)

where 5730 is the ¹⁴C half-life, A is the ¹⁴C concentration (as percent of modern carbon) in the groundwater, and A₀ is the fraction of the current total carbon concentration, attributable to carbon in the recharge area.

Table 1. ¹⁴C-based ages of groundwater from selected wells in the Lower Chao Phraya (LCP) Basin [45].

Well No.	Model Layer No.	Aquifers Name	Observed Groundwater Age (Years)
1	3	PD	7100
2			9100
3			7300
4			18,100
5			5300
6			6900
7			9500
8			9800
9			5200
10			8300
11			18,300
12			6800
13			10,000
14			7700
15			8500
16			5700
17			2300
18	4	NL	8100
19			7400
20			8900
21			18,900
22			18,100
23			5500
24			20,900
25			19,400
26			19,900
27			12,500
28			8600
29			8400
30	5	NB	6300
31			9000
32			8500
33			18,200
34			1600
35			13,800
36			9100
37			154,000

shows the ¹⁴C-based ages of groundwater from selected wells in the LCP Basin [45]. The groundwater age was 8500 (PD), 13,000 (NL) and 10,200 (NB) years, with mean ~11,000 years. This showed that groundwater is “young” (>50 years and <12,000 years) or fossil groundwater (>12,000 years), following the labels applied by Bierkens and Wada, Gleeson et al., and Jasechko et al. [50,64,65].

3. Methods

3.1. Model Concept and Boundary Conditions

The model used was described by Tanachaichoksirikun et al. [60] and derived from a database that included hydraulic parameters (hydraulic conductivity and specific storage), groundwater budgets (recharge, pumping and river), topographic data (Digital Elevation Model: DEM and aquifer thickness) and boundary conditions (no-flow boundary and constant head), that was calibrated by matching groundwater hydraulic head, hydraulic gradient and sensitivity of recharge and pumping in 2009–2014. The parameters used were the same as those used for verification in 2007–2008.

Figure 3 shows a 3D groundwater model of the LCP basin. The groundwater model was separated into nine layers—eight for groundwater aquifers and a topmost layer for Bangkok Clay (BKC). A constant head was assigned in the north side for groundwater flow from the Upper Chao Phraya basin and in the south for the tidal water of the Gulf of Thailand, using the hydraulic head near the edge. River boundaries were applied for four rivers using data from 27 river monitoring stations and calculated following Harbaugh et al. [52]. The recharge boundaries were defined by land use, soil property and rainfall intensity [4,23]. Below the lowest layer, there was a no-flow boundary of impermeable bed rock. The west edge had a no-flow boundary for the Tenasserim hills, and the east had a no-flow boundary defined by small hills. In addition, the major groundwater pumping in PD, NL, NB aquifers—due to their heavy consume groundwater for domestic, industry and agriculture and many available data for calibration and validation—were treated as pumping. For these reasons, this research focused on the groundwater age distribution in these aquifers.

3.2. Hydraulic Parameters

Hydraulic parameter data were derived from our previous study [60]: we considered a hydraulic gradient technique to estimate hydraulic parameters, hydraulic conductivity and anisotropy and concluded that the anisotropy of the LCP basin was 10⁴. In addition, here, we considered the porosity using the observed groundwater age from data of Sanford and Buapeng [45] and followed Meyer et al.’s technique [66]. The groundwater age data were used for verifying porosity of the LCP basin at the sampled points. Backward tracking was applied from the derived age in the observation well to the source, to compare with the model prediction. The groundwater parameters are shown in

Table 2. Groundwater parameter in the LCP basin.

Aquifers	Porosity	Kh (m/s)	Kv (m/s)	Ss (m−1)
Bangkok Clay and Unconfined	0.03	1 × 10−8–1 × 10−9	1 × 10−8–1 × 10−9	0.03–0.35 (Sy)
Bangkok	0.2–0.3	5 ×10−5–6 × 10−6	5 × 10−9–6 × 10−10	1.10 × 10−5–4.80 × 10−5
Phra Pradeang	0.25–0.35	2 × 10−5–7 × 10−6	2 × 10−9–7 × 10−10	3.35 × 10−5–4.50 × 10−6
Nakorn Luang	0.2–0.35	1 × 10−4–6 × 10−6	1 × 10−8–6 × 10−10	3.85 × 10−5–9.50 × 10−6
Nonthaburi	0.3–0.35	2 × 10−4–7 × 10−6	2 × 10−8–7 × 10−10	3.95 × 10−5–9.50 × 10−6

.

3.3. Future Climate Scenarios

The Intergovernmental Panel on Climate Change (IPCC) [8] developed a new greenhouse gas model, the Representative Concentration Pathways (RCP), reported in The Fifth Assessment Report (AR5). RCP was developed by using an initial greenhouse gas concentration and predicting the impact on climate change for several greenhouse levels. RCPs used here were based on four scenarios, i.e., RCP2.6, RCP4.5, RCP6.0 and RCP8.5. The numbers, 2.6, 4.5, 6.0 and 8.5, are labels for a set of representative pathways greenhouse gas concentrations defined in the IPCC AR5 report.

We assumed that climate change was reflected in changes in rainfall and that the recharge rate was directly related to rainfall. The projected climate scenario during 2020–2099 was derived from Wattanasetpong et al. [67], based on records for 30 years of Thailand’s LCP monitoring stations and the Royal Netherlands Meteorological Institute (KNMI) Climate Explorer [68]. They considered the six important criteria for climate change—rainfall, temperature, pressure, humidity, evaporation and water discharge to the ocean. Ruangrassamee et al. [69] investigated and showed that the IPSL-CM5A-MR (Institute Pierre Simon Laplace) scenario was optimal for the Chao Phraya water shed, because the climate model was a minimum root-mean-square error and no bias data. Additionally, three Representative Concentration Pathways (RCP) values were selected—2.6, 4.5 and 8.5—as representatives of a range of carbon dioxide concentrations and rainfall rates, where RCP 4.5 represents rainfall close to the present [67]. RCP 6.0 was not considered, because it represents an intermediate scenario between RCP4.5 and RCP8.5.

Figure 4 illustrates the predicted future climate from downscaling projected annual rainfall from IPSL-CM5A-MR under the RCP 2.6, 4.5 and 8.5 models, from 2020 to 2099, compared to a baseline trend [67]. The baseline period was based on rainfall from 1984 to 2014 for checking the change in groundwater age distribution with climate change. In the climate scenarios, monthly precipitation was predicted for 2020 to 2099 [67]. We predicted groundwater sustainability in the basin by analyzing movement of groundwater from the groundwater age distribution. Note, future groundwater extraction assumed a pumping rate that was the same as the baseline, because the area showed a very slowly increasing pumping rate.

The annual rainfall trends predicted by the various model are shown in Figure 4. The annual rainfall of the RCP 2.6 model showed a gradually increasing trend, whereas the RCP 8.5 annual rainfall had a significantly increasing trend from 2020 to 2099. However, the RCP 4.5 model showed a more slowly increasing trend than the RCP 8.5 model, in the period 2036–2099. IPSL-CM5A-MR with RCP 8.5 predicted the highest rainfall (↑7.4%), whereas IPSL-CM5A-MR RCP 2.6 predicted the smallest (↓17.4%). Close to the current trends were found in IPSL-CM5A-MR with RCP 4.5 (↑3.1%).

3.4. Predicted Groundwater Age

We focused on change in groundwater age under global climate change. The particles, assumed to be groundwater samples, were randomly assigned to sampling points in the PD, NL and NB aquifers. The predicted groundwater age distributions were simulated, in every climate scenario using the particle backward tracking method, because we had known sample locations and could track backwards from them. However, we analyzed on a regional scale, because groundwater came from several sources, e.g., mixing between young and old groundwater, pumping and recharge, so the particle was used to infer the geographic distribution of groundwater ages.

For model evaluation, the dataset from 2007–2014 was divided into a calibration period (2009–2014) and a validation period (2007–2008). The MODFLOW-2000 model was calibrated from MODFLOW’s parameter estimation and also manual parameter estimation. The model performance can be considered satisfactory if absolute residual mean error ≤ 5 m, NRSE ≤ 10% and R² > 0.85 [70]. For model prediction, the 2020–2099 data were projected to predict the groundwater age distribution, using initial particles in the calibration model.

4. Results and Discussion

4.1. Calibration and Verification

The calibration and verification by Tanachaichoksirikun et al. [60] for the correlation of simulated and measured hydraulic head is shown in Figure 5. In Figure 5a, the simulated and measured hydraulic heads were compared, in 2009–2014, with an absolute residual mean error of 4.6 m and NRSE = 5.9% and R² = 0.91.

In the verification model (Figure 5b), the simulated and measured hydraulic heads were compared in 2007–2008. The absolute residual mean error was 2.6 m, with NRSE = 3.8% with R² = 0.95. So, our errors were less than those of Du et al. [70], often used as a benchmark for numerical simulation of groundwater flow: for any node, the absolute residual mean error should be less than 5 m, with NRSE < 10%.

In addition, here, we calibrated using change in hydraulic gradient to anisotropy as in our previous work [60]. The calibration using only the hydraulic head was not sufficient in a regional scale, because the LCP basin is alluvial. The aquifers are sedimentary units, with high anisotropy. the data showed that the optimal anisotropy was 10⁴, due to the high correlation between hydraulic head and hydraulic gradient (R² > 0.8). Additionally, this model was analyzed for sensitivity to recharge and pumping and showed that groundwater flow was comparable with the real-world condition.

4.2. Simulated Versus Observed Groundwater Age

Figure 6 shows a close agreement between the observed ¹⁴C groundwater ages from Sanford and Buapeng [45] (over 11,000 years) and our simulation (over 10,733 years). The standard deviation was 3300 years, with R² = 0.81. Thus, this model can be used for predicting groundwater age distribution: R² > 0.8 and standard deviations in the range of 2000–4000 years [45].

4.3. Predicted Groundwater Age Distribution

Figure 7 illustrates the groundwater age distribution in the LCP basin, for the base case and IPSL-CM5A-MR with RCP 2.6, 4.5 and 8.5. The groundwater age was shown separately for each aquifer and overall. In general, the younger groundwater was near the recharge and injection zones, where the groundwater age is zero and increased with gradient and depth. The groundwater age in the LCP basin was distributed widely, i.e., ages varied over a wide range; because of the regional basin, some groundwater was resident for a long time, travelled long distances and several sources of water were involved—as discussed in the previous section.

Under the base case (Figure 7a), groundwater age in the PD aquifer covered 100 to 10,000 years, and the mean groundwater age was ~2800 years. The groundwater age in the NL aquifer covered 100–70,000 years, and the mean age was ~8000 years. The NB aquifer distribution covered 1000–100,000 years with a mean age of ~14,500 years. Overall, the average age was ~11,000 years. Thus, the study area was mainly old groundwater. Groundwater flowed slowly, because the basin was a sedimentary aquifer and clay layers between the aquifers hindered groundwater flow from the surface area to the aquifer.

Under RCP2.6 (Figure 7b), groundwater age in the PD aquifer was 100 to 9000 years, with average age = 2988 years (i.e., 188 years longer than the base case). The average age distribution was 67% wider than the base case. Groundwater age in the NL aquifer was 200–80,000 years, with an average of 8986 years (+985 years from the base case). The age distribution was also wider, by 80.1%, compared to the base. Groundwater age in the NB aquifer was 1000–120,000 years. The average age = 18,253 years (+3745 years from the base). The age distribution was 62% wider than the base case. This showed that the lower rainfall contributed to increased groundwater age. The lower rainfall contributed to decreased recharge, low groundwater budget, which, in turn, led to lower hydraulic head and lower velocity. The shallow aquifer was affected by the climate change, due to lower distance and thinner clay layers.

Under RCP4.5 (Figure 7c), groundwater age in the PD aquifer was 100–8000 years with average = 2850 years (+50 years). The distribution was 59% wider. Groundwater age in the NL aquifer ranged from 100 to 70,000 years. The average age was 7967 years (−34 years, i.e., less than the base case). The distribution was wider by 83.2%. In the NB aquifer, groundwater age covered 800–120,000 years with average = 14,911 years (+411 years). The distribution was 59% wider. This showed that the predicted rainfall was close to the base case. The groundwater age was similar to the base case.

Under RCP8.5 (Figure 7d), groundwater age in the PD aquifer covered 100–8000 years with the average = 2746 years (−69 years). The distribution was 80% wider. Groundwater age in the NL aquifer covered 100–60,000 years. The average age = 8914 years (+913 years). The distribution was 89% wider. In the NB aquifer, groundwater age covered 1000–120,000 years, with average = 17,378 years (+2870 years). The distribution was 62% wider. This showed that the higher rainfall led to decreased groundwater age in a shallower aquifer, but in a deeper aquifer, groundwater age was apparently still increasing, but our simulation time extended only to 100 years. The rainfall had not flowed to the NB aquifer, so the groundwater age in NB was increased.

Additionally, the average groundwater age in a shallow aquifer was younger than in a deeper one. The groundwater age distribution in the PD aquifer covered 100–10,000 years, whereas, in the NL and NB aquifers, it extended to 100,000 years.

Although Engdahl and Maxwell [49] suggested that climate change affected groundwater age, our work showed that groundwater age changed only slightly in transient groundwater flow conditions, because the climate change was not accompanied by high rainfall, at that time. In higher recharge scenarios, e.g., IPSL 8.5, the groundwater age distribution shifted slightly, because rainfall increased groundwater velocity and rain generated more surface water. In contrast, in lower recharge scenarios, e.g., IPSL 2.6, the groundwater age shifted significantly, because lower rainfall caused longer residence times and reduced flux (slower velocities).

Overall, changing groundwater age distributions, with differing climate change scenarios, reflected many changes of hydraulic head, recharge and velocity. The change in the age distributions showed the dynamics of the groundwater systems, as fractions of groundwater shifted slightly at various ages, from 10 to 100 years as shown in Figure 7. The age distribution changes indicated changes in location and rate of recharge. Then, groundwater levels, that cause land subsidence, may be unstable. In climate change scenarios, with increased rainfall, making the regional system wetter, the upper aquifer was more completely filled and the excessive groundwater flowed to the surface water system. Groundwater age distributions were unaffected. Change in water budget between surface and groundwater was reduced and caused an increase in older groundwater, although the groundwater ages rarely changed.

4.4. Distribution of Groundwater Age under Climate Change Impact

The projected climate change, under the IPSL-CM5A-MR, was simulated to check for future groundwater age. More recharge, therefore, increasing groundwater levels, led to younger groundwater in shallow aquifers and older water in deeper aquifers. Thus, future groundwater recharge is the most significant feature that controls decrease or increase of age in each aquifer. Increased future groundwater age and the time to recharge was longer than human life spans, which indicated that groundwater in the LCP basin is nonrenewable groundwater. The groundwater recharge affected the groundwater fluctuation, which led to the groundwater problems, e.g., land subsidence and seawater intrusion.

In the shallower aquifer, groundwater ages decreased because, as older groundwater was pumped out of the lower aquifer, it was replaced by younger groundwater from the upper layer. The old groundwater that was pumped out of the aquifer cannot be replaced by recharge, because the effects of old groundwater were replaced by young water in the upper part of the aquifer and mixed with low velocity.

Nonrenewable groundwater will contribute to land subsidence and consequent detrimental effects, e.g., damage to infrastructure and increased flood damage and risk due to lower ground surface level [71]. Soil is an elastic material, so if groundwater is recovered, the soil level will be increased. However, if groundwater level decreased for long periods and contributed to land subsidence for longer times, soil material will creep and become inelastic, so although groundwater level is restored, the land will not return to previous levels.

Since the LCP basin connects to the Gulf of Thailand, pumping groundwater has led to a progressive increase in salinity: decrease in hydraulic head, under the influence of pumping, especially in low recharge scenarios, led to rise of the sea–groundwater interface, leading to more saline groundwater resources. This effect was more obvious in the shallower aquifers, because of the higher sensitivity to hydraulic heads and to higher pumping, induced by lower costs.

The predicted increase in groundwater age from recharge will affect the groundwater level. Therefore, groundwater pumping must be carefully managed, because it may cause land subsidence, seawater intrusion and salinity. Moreover, pumping in nonrenewable groundwater may contribute to groundwater shortage, because rainfall cannot recharge it in human life spans. Here, we showed that the groundwater age distribution, in the deeper layers, i.e., the NL and NB aquifers, will gradually increase in the order of 100 years, indicating that recharge will not fill these aquifers; therefore, pumping in these aquifers must be carefully managed.

5. Conclusions

Variations in precipitation and recharge affect groundwater age distributions on regional scales. We assessed the distributions of groundwater age, under climate change in several scenarios, in Thailand’s Lower Chao Phraya (LCP) basin. The optimal climate scenario, IPSL-CM5A-MR, derived in work by Ruangrassamee et al. and Wattanasetpong et al. [67,69], was used here, and we varied the Representative Concentration Pathway (RCP) models (labeled 2.6, 4.5 and 8.5), from 2020 to 2099. Simulations used MODFLOW-2000 and MODPATH under transient-state conditions. We found that groundwater age slightly decreased in the shallower aquifers, whereas the age increased in the deeper aquifers, under increasing recharge scenarios, because current rainfall did not move water to the deep aquifers. The overall average groundwater age gradually increased, due to groundwater age mixing in both shallow and deep aquifers. However, under decreasing recharge scenarios, groundwater age increased in both shallow and deep aquifers, because reduction of recharge caused groundwater to flow slower and for longer times, so that groundwater became older, which was also reflected in decreased in hydraulic heads.

The fluctuation of groundwater recharge led to change in groundwater age, which was an indicator of the groundwater level, which, in turn, allowed for land subsidence, seawater intrusion and increased salinity. Further, groundwater became short in the area, if groundwater was over pumped, especially in the NL and NB aquifers, where groundwater age increased in both increased and decreased recharge scenarios.