The Inﬂuence of Climate Variability Effects on Groundwater Time Series in the Lower Central Plains of Thailand

: This research studies the relationship between the climate index and the groundwater level of the lower Chao Phraya basin, in order to forecast the groundwater level in the studied area by using Autoregressive Integrated Moving Average with Explanatory (ARIMAX). The combination of 6 climate indices—Dipole Mode Index, Indian Summer Monsoon Index, Multivariate ENSO Index, Sea Surface Temperature NINO4, Southern Oscillation Index and the Western North Pacific Monsoon Index—were used, along with the groundwater level data from 14 stations during the period 1980–2011 to develop the forecast model and verify it with the data of 2012.The first step was correlation of the ARIMA model with Autocorrelation Function and Partial Autocorrelation Function. The possible model was then selected using BIC statistics. Diagnostic Checking was done to consider the white noise characteristic of estimated residuals by using the statistics of Box and Ljung (Q-statistic). If the selected models were found to be proper, then the Granger Causality Test of the leading parameters or the climate index would be performed as the next step. The results show that there is a relationship between the groundwater level and the climate index. The model could be used to forecast effectively the average RMSE value at 0.6. The last step was to develop the MODFLOW for a conceptual model and synthesize groundwater levels in the study area, which covers around 43,000 km 2 and has 8 layers of groundwater, with Bangkok clay on the top. All other boundary values were set to be steady. The calibration was done using the data of 325 observed wells. The normalized RMS value was 9.705%. The results were verified by the data using ARIMAX over the same time periods. To conclude, the simulated results of the monthly groundwater level in 2012 of the wells have a confidence interval of around 95%, which is near the result from the ARIMAX model. The advantages of the ARIMAX model include high accuracy, no requirement for a large amount of data and inexpensive implementation. It is one of the effective tools for the groundwater prediction.


Introduction
As is widely known, the world's weather is currently facing global changes [1], which can be seen from the variability of weather and changes in the climate in many parts of the world [2]. Every aspect is involved in this climate variability, including atmosphere, ocean levels and average temperature, resulting in changes in season timing and rain distribution [3]. The studies of the Intergovernmental Panel on Climate Change (IPCC) in 2001 [4] showed that the average world temperature had increased by 1.4 to 5.8 degrees Celsius. This resulted in the melting of polar ice and the expansion of the ocean area, the increase of sea level by at least 0.09 m to 0.88 m and variability in the distribution of rain, which could result in both more severe drought and flooding.

Linkages between Climate Variability and Water Resources
Worldwide, people's water resources respond to climate variability very quickly. Regional climates are affected by many forms of climate variability. The one that is most well-known is the El Nino-Southern Oscillation (ENSO) [7]. Moreover, variations in the phenomena were found in decade-long spans such as Pacific Decadal Oscillation (PDO) [8,9], the North Atlantic Oscillation (NAO) [10] and the Atlantic Multidecadal Oscillation (AMO) [11].
Previous studies have shown the relationship between rainfall and climate variability in many areas throughout the world, such as in Queensland by Chiew et al. [12] and in the southwestern region of the United States by Hanson et al. [2]. Thailand has been affected by the difference in climate from both the Indian and Pacific Ocean. Yearly rainfall decreases in El Nino years and increases in the La Nina years [13]. The relationship between climate index and change of multi-year rainfall, up to decades, was researched and found by Kusreesakul [14], Limsakul and Goes [15]. The variations in the world's climate have also impacted surface water resources and groundwater resources in Thailand. Limsakul et al. [16] has found that the yearly raining period is related to the Indian Ocean Monsoon Index (IMI) and the Western North Pacific Monsoon Index (WNPMI).
In the future, water use will rely on groundwater more because of the quickly changing and unreliable surface water level. In many areas, the surface water level is forecast to become more variable, which could drop the quality of surface water because of more severe drought and more frequent rain [17]. In the IPCC report, current surface water management is not sufficient to ensure the quantity of tap water impacted by climate change. Therefore, the number of wells is increasing very quickly [18] because farmers rely on it more and the surface water supply will be unreliable due to climate change.
Although the number of journal publications per year and the cumulative number of papers about the impact of climate change on ground water resources from 1990 to 2010 had largely increased [1], more research is still needed urgently to define conjunctive use adaptation strategies to improve the management of groundwater [19], under the impact of world climate change and variation [20]. Studying the probable impact of climate change on groundwater variation is much more complex than the impact on the surface water [21]. The groundwater could have a lifetime from a few days to a hundred thousand years. Although groundwater could slow down the impact from climate change, the impact on the groundwater is difficult to monitor and detect [22]. Moreover, many activities caused by humans, such as pumping up the groundwater, could affect the groundwater in the same period as the climate change, which makes the differentiation of impact from climate change and human activity more complex. USGS has tried to study the groundwater level hydrology and agronomical chemistry response signal in the interannual to multidecadal time scales because this time scale variation could be highly significant for groundwater resource management [2,[23][24][25].
The main interest of groundwater and the climate change research is the forecasting and estimation of probable direct impact on the rainfall variation and temperature formation [26][27][28]. These studies could be simulated by different models, such as the water and soil equilibrium model [28,29], Empirical model [30], conceptual model [31] and complex distribution model [32,33]. In addition, Changnon et al. [34] have studied the relationship between monthly precipitation (P) and the shallow groundwater level (GW) in 20 wells scattered across the American state of Illinois using autoregressive integrated moving average (ARIMA) modeling. The model found that a lag of 1 month between P to GW has the strongest temporal relationship found across Illinois, followed by no (0) lag in the northern two-thirds of Illinois, where mollisols predominate and a lag of 2 months in the alfisols of southern Illinois. Adhikary et al. [35] also used ARIMA to simulate the fluctuation of the groundwater table of all monitored wells in the Kushtia district of Bangladesh. The results show that the predicted data represented the actual data very well for each monitored well. Furthermore, groundwater heads at a confined aquifer in southwest Florida show a nonstationary long-term (multi-year) fluctuation. Ahn and Salas [36] introduced an approach to build a time series models of nonstationary data at different time intervals based on an observed time series sampled at a reference interval. The model utilized in their study was a first-order difference ARIMA model. However, some groundwater head data may also be fitted adequately by a second-order difference time series model [37]. Five-time series models were applied: autoregressive (AR), moving average (MA), auto-regressive moving-average (ARMA), autoregressive integrated moving-average (ARIMA) and seasonal auto-regressive integrated moving-average (SARIMA). The results showed that the AR model with a two-time lag (AR(2)) shows the best forecasting of groundwater level. However, they should be combined with several time series models for a better prediction of groundwater level [38].
The MODFLOW model, which was developed by the U.S. Geological Survey, is widely used for groundwater modelling [39]. Visual MODFLOW was then developed by Waterloo Hydrogeologic, Inc. (Kitchener, ON, Canada) to help prepare the data and conveniently display the groundwater model. For instance, Lawrence et al. [40] used MODFLOW for the hydrogeology study in the Hat Yai district of Songkhla province in Thailand to see the influence on the amount of groundwater usage in an urban area. Margane et al. [41], with the Department of Mineral Resources, applied MODFLOW in the hydrogeology study in the Chiang Mai-Lampoon basin to organize the quality of groundwater data and to show the risk of contaminated groundwater. The Department of Groundwater Resources (DGR) [42] did research about the addition of water to the groundwater through the pond system in the northern watershed area (Phitsanulok, Sukhothai and Phichit) and they did a study on the impact of underground structure due to the restoration of water pressure in Bangkok and vicinity. Moreover, choosing the parameters in different modeling programs, such as hydraulic conductivity, is an important step [43,44], because any groundwater basin might have many k values, since the soil and rock layers are different. Most of the modeling on Thailand data used the anisotropic ratio at about 10 times [44]. In a study similar to this research study, the lower Chao Phraya River basin and the surrounding area was also researched by Arlai et al. [45], who looked into the effective approaches to assessing the reliable parameters in the Bangkok aquifer model. The previous research is greatly beneficial to this study. However, using a lot of surrounding area was also researched by Arlai et al. [45], who looked into the effective approaches to assessing the reliable parameters in the Bangkok aquifer model. The previous research is greatly beneficial to this study. However, using a lot of data may take a long time to do modeling and the scope of the model may be limited. The data assigned for the model ought to be the minimum necessary, for the sake of simplicity [46].
Climate change impacts the groundwater resources critically. Other change and trends that are difficult to forecast also have an important impact that cannot be neglected. In particular, climate change can impact the quality and quantity of groundwater more and more, so prediction models are important, because groundwater could become a more important and reliable water resource than surface water.

Study Area
The groundwater basin in the lower Chao Phraya River is located in one of the most important and populous parts of Thailand, the lower central region. As shown in Figure 1, the area covers about 20 provinces with approximately 43,333 km 2 and has the biggest groundwater basin in Thailand, at 269.312 billion m 3 . The area can be divided into 2 major parts. The lower plains part covers the area from the plain in Manorom city in Chainat province to the estuary area of the Chao Phraya River, which contains the water from the ground level to about 600 m deep. The other part is the borders of the Chao Phraya basin on both the east and west side, which is a mountainous area and has less water. The west side of the Chao Phraya basin edge covers the area from Uthai Thani province and the west side of Suphanburi province to Nakorn Pathom province. The east side of the Chao Phraya basin edge covers the area from Lopburi province, Saraburi province, Nakorn Nayok province, Prachinburi province and Cha Choeng Sao province. Most of the area is the lower plains, except the north of Nakorn Pathom province, which is an upland area. The Chao Phraya River is the major river that runs from the north of Phra Nakorn Sri Ayutthaya province to the estuary area in Samut Prakarn province. The characteristics and types of Most of the area is the lower plains, except the north of Nakorn Pathom province, which is an upland area. The Chao Phraya River is the major river that runs from the north of Phra Nakorn Sri Ayutthaya province to the estuary area in Samut Prakarn province. The characteristics and types of rock and the geological structure of the edge of lower central region plain are mostly river sediment that flowed into the sediment basin, which is composed of Quaternary deposits. The hydrogeology of the lower central region plains is composed of both unconsolidated rocks and consolidated rocks [47].
From the drill hole in the studied area, which collected data to a depth of 700 m, it was found that there is a groundwater level, which is the gravel level and which could be divided into 8 levels, as shown in Figure 2, which is modified from Asian Institute of Technology [48]. rock and the geological structure of the edge of lower central region plain are mostly river sediment that flowed into the sediment basin, which is composed of Quaternary deposits. The hydrogeology of the lower central region plains is composed of both unconsolidated rocks and consolidated rocks [47].
From the drill hole in the studied area, which collected data to a depth of 700 m, it was found that there is a groundwater level, which is the gravel level and which could be divided into 8 levels, as shown in Figure 2, which is modified from Asian Institute of Technology [48]. The water levels in wells that are screened in the upper-most unit (Bangkok aquifer) are going to respond differently from the water levels in wells in deeper hydrogeologic units. Because the upper-most unit wells have been influenced by the ground water recharge such as precipitation and stream runoff, the response happens more quickly than in other wells.

Methodology
The data collected included a climatic/oceanography index and groundwater levels, to study the interrelation between them by analysis of ARIMAX. The water level forecasting of the ARIMAX method was compared with groundwater levels obtained from the simulation MODFLOW model, a flowchart of methodology shown in Figure 3 and the details are discussed below. The water levels in wells that are screened in the upper-most unit (Bangkok aquifer) are going to respond differently from the water levels in wells in deeper hydrogeologic units. Because the upper-most unit wells have been influenced by the ground water recharge such as precipitation and stream runoff, the response happens more quickly than in other wells.

Methodology
The data collected included a climatic/oceanography index and groundwater levels, to study the interrelation between them by analysis of ARIMAX. The water level forecasting of the ARIMAX method was compared with groundwater levels obtained from the simulation MODFLOW model, a flowchart of methodology shown in Figure 3 and the details are discussed below.

Climatic Indices Data
According to NOAA, the major climate variation of this region is the Asian summer monsoon, which is analyzed with the Indian Summer Monsoon Index (IMI) and Western North Pacific Monsoon Index (WNPMI) and the Indian Ocean Dipole (IOD) is analyzed with the Dipole Mode Index. Those indices were used to consider the differences in sea water surface temperature around the equator area from the western Indian Ocean (50-70 • E and 10 • S-10 • N) to the southeastern part of the Indian Ocean (90-110 • E and 10 • S-10 • N).
Another phenomenon in this region is the El Nino-Southern Oscillation (ENSO), which involves several indices, such as the Multivariate ENSO Index (MEI), Southern Oscillation Index (SOI) and Sea Surface Temperature (SST). NINO4 is in the western equatorial part of the Pacific Ocean, or the so-called warm water basin and is located in latitude 5 • N-5 • S and longitude 150 • W-160 • E, which is the area that has the highest sea water temperature. It can affect many countries in the western Pacific; a total of 6 indices are used in analysis. The data also has been collected from 1960 to the present.

Groundwater Data
The monthly groundwater level data was monitored in wells in the lower Chao Phraya basin by the DGR. The data was collected from 14 wells in the studied area, which includes about 15 provinces; information from 4 sample stations is shown in Figure 4.  From the groundwater data that was collected from 1980 to 2012, measured from the land surface at the CT48-2 station, the water level has the highest variation and at the CT23 station, the water level has the lowest variation. At the same time, the average water level at the edge of the well is highest at CT33-2 station, at 47.60 m. The overall groundwater level has the negative kurtosis and negative skewness distribution.

Theory
When all necessary data, such as climate/oceanography index and groundwater level, had been collected, the model was developed using ARIMAX. It is one of the most popular statistical models in economic analytics for many types of time series forecasts, such as monthly incoming visitor forecasts [49] and daily car count forecasts [50]. ARIMAX has been constantly developed from the ARIMA mode, which had originally been presented by Box and Jenkins [51].
The climate/oceanography index+ and groundwater level used in the analysis should be stationary, so that the time series data has the constant mean (µ), variance (σ 2 ) and covariance (γ). Therefore, the characteristics of the data need to be inspected with the unit root test, which can be done by using the Augmented Dicky-Fuller Test (ADF) [52] in 3 forms, as shown in Equations (1)-  From the groundwater data that was collected from 1980 to 2012, measured from the land surface at the CT48-2 station, the water level has the highest variation and at the CT23 station, the water level has the lowest variation. At the same time, the average water level at the edge of the well is highest at CT33-2 station, at 47.60 m. The overall groundwater level has the negative kurtosis and negative skewness distribution.

Theory
When all necessary data, such as climate/oceanography index and groundwater level, had been collected, the model was developed using ARIMAX. It is one of the most popular statistical models in economic analytics for many types of time series forecasts, such as monthly incoming visitor forecasts [49] and daily car count forecasts [50]. ARIMAX has been constantly developed from the ARIMA mode, which had originally been presented by Box and Jenkins [51].
The climate/oceanography index+ and groundwater level used in the analysis should be stationary, so that the time series data has the constant mean (µ), variance (σ 2 ) and covariance (γ). Therefore, the characteristics of the data need to be inspected with the unit root test, which can be done by using the Augmented Dicky-Fuller Test (ADF) [52] in 3 forms, as shown in Equations (1)- (3).
The parameter that is considered in every equation is δ. If δ = 0, Y t will have a unit root from comparing the calculated t-statistic value and the value from the Dickey-Fuller tables, or the MacKinnon critical values.

ARIMA Model
The Box-Jenkins method is the way to find the proper model for the time series value by using the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) to consider the model. The model is composed of 3 major parts, which are the Autoregressive Model (AR(p)), Integrated Process (I(d)) and Moving Average Model (MA(q)).
Autoregressive (AR(p)) is the model to show the observed y t value; it was assigned from the values of y t−1 ,...,y t−p or the previous observed value (p) by using the process of system AR(p), which is the process of the Autoregressive system at p order, which can be written as shown in Equation (4).
where y t is the groundwater level at time t; µ is the constant value; φ i is the I ordered parameter; and ε t is the error at time t. Moving Average (MA(q)) is the model to calculate the observed y t value, which is assigned from the error at ε t−1 , . . . ,ε t−q or the previous error by using the process or system MA(q), which is the Moving Average at the q order, which can be written as in Equation (5).
The last step is to assign the proper model of the ARIMA model, which can be done by considering the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF). In considering the Autoregressive Process model at p ordered and Moving Average at q ordered would result in Autoregressive Integrated Moving Average mode at p, d and q ordered, or ARIMA (p,d,q), which can be written as shown in Equation (6).
where y t is the groundwater level at time t; p is the order of Autoregressive; d is the number of times that differences were found in the constant characteristic time series data; q is the order of Moving Average; ∆ d is the sign of finding the difference at d ordered (d-th difference operator); δ is the constant value; φ 1 , φ p is the parameter of Autoregressive; θ 1 , θ q is the parameter of Moving Average; ε t is the error at the time t, which has been assigned to be white noise; and t is the time index.

ARIMAX Model
ARIMAX is the model that was adapted by including X, which is the leading variable (or Structural Variable or Exogenous Variable) to improve the forecasting accuracy. In this study, this variable is the climate and oceanography index and was newly formulated into Equation (6) to create Equation (7), as follows.
When the independent variable is included in the ARIMA model, it results in the ARIMAX model, shown in Equation (8).
In Equation (8), x t is the climate/oceanography index, where θ(L) is the degree of polynomial of the Exogenous variable that affects the Endogenous variable.
In testing to see whether X influences Y or not but using the null hypothesis (that is, H 0 : X does not influence Y) which can be said that to conclude that X influences Y, the null hypothesis needs to be rejected and when the variable passes the leading index testing, that variable could then be used to analyze the ARIMAX model.
When the proper order of ARIMAX model has been set up from the consideration of Bayesian Information Criteria (BIC), diagnostics should be done by the insignificant Ljung-Box Q statistic.
After the most proper model is devised, then the groundwater level should be forecast to compare with the observed value, to compute the Root Mean Square Error (RMSE) and Mean Absolute Percent Error (MAPE). Finally, the groundwater level forecast by the statistical process would then be compared to the groundwater level that was synthesized from the physical model, which is MODFLOW in this case.

Conceptual Groundwater Model
The analysis of data composed of geological data, hydrogeological data, groundwater characteristics, geographical data, groundwater usage data, mapping data and hydrological cross section graphics, which are partly modified from Arlai et al. [45], Department of Groundwater Resources [47] and Fornés and Pirarai [53], resulted in the conceptual model by MODFLOW as shown in Figure 5. The model is rectangular between (north-south) UTM 1,450,000 to 1,800,000 m and (east-west) 500,000 to 800,000 m. By having a fixed grid size at 1 km 2 , the model has a total of 945,000 grids, which are divided into 9 layers, comprising 8 layers of groundwater and the average 20-m thick Bangkok clay as the top layer.
Water 2018, 10, x; doi: FOR PEER REVIEW 9 of 24 compared to the groundwater level that was synthesized from the physical model, which is MODFLOW in this case.

Conceptual Groundwater Model
The analysis of data composed of geological data, hydrogeological data, groundwater characteristics, geographical data, groundwater usage data, mapping data and hydrological cross section graphics, which are partly modified from Arlai et al. [45], Department of Groundwater Resources [47] and Fornés and Pirarai [53], resulted in the conceptual model by MODFLOW as shown in Figure 5. The model is rectangular between (north-south) UTM 1,450,000 to 1,800,000 m and (eastwest) 500,000 to 800,000 m. By having a fixed grid size at 1 km 2 , the model has a total of 945,000 grids, which are divided into 9 layers, comprising 8 layers of groundwater and the average 20-m thick Bangkok clay as the top layer.

Boundary
The hydrogeological model employs the following assumptions. The groundwater in the study area flows from north to south, which shows that horizontal permeability is better than vertical. The density of water is constant. Each aquifers layer is heterogeneous and anisotropic. The confining beds are homogeneous and isotropic. The groundwater flow is in steady state and the no-flux boundary around the area means that all aquifers are inactive cells, as shown in Figure 6, except for the south

Boundary
The hydrogeological model employs the following assumptions. The groundwater in the study area flows from north to south, which shows that horizontal permeability is better than vertical. The density of water is constant. Each aquifers layer is heterogeneous and anisotropic. The confining beds are homogeneous and isotropic. The groundwater flow is in steady state and the no-flux boundary around the area means that all aquifers are inactive cells, as shown in Figure 6, except for the south set, which is the constant head boundary. The water level at the Gulf of Thailand edge is set to be +0.00 m (MSL), according to Barlow [54] and Langevin [55], as shown in Table 1. The proper boundary of the model is in accordance with the hydrogeological characteristic of groundwater basin system in lower Chao Phraya from the study of DGR, such as the aquifer type of the Bangkok clay, on the top, which is set as an unconfined aquifer. But the groundwater from the Bangkok (BK) aquifer all the way to the Pak Nam (PN) aquifer is set as confined.

Parameters
The Digital Elevation Model (DEM) data is imported in order to have the height of the top and bottom of each layer as the model. The ground surface elevation is imported from the elevation data of the DGR national chart at the scale of 1:50,000. The study area has the elevation level at around 0-900 m above MSL. The initial head of every water layer is set from the groundwater level using the data from the wells observed by the DGR in 2009.
The assumption initially set the soil layer to be heterogeneous and anisotropic, referred from the hydraulic conductivity in the model by Spitz and Moreno [56]. The initial values for the calculation are shown in Table 2. Table 2. Initial hydraulic conductivity, K, for model.

Parameters
The Digital Elevation Model (DEM) data is imported in order to have the height of the top and bottom of each layer as the model. The ground surface elevation is imported from the elevation data of the DGR national chart at the scale of 1:50,000. The study area has the elevation level at around 0-900 m above MSL. The initial head of every water layer is set from the groundwater level using the data from the wells observed by the DGR in 2009.
The assumption initially set the soil layer to be heterogeneous and anisotropic, referred from the hydraulic conductivity in the model by Spitz and Moreno [56]. The initial values for the calculation are shown in Table 2. Table 2. Initial hydraulic conductivity, K, for model.

River
River Stage (m. msl.) Riverbed Bottom (m. msl.) The data of the major rivers in the study area-the Chao Phraya River, Tha Chin River, Pa Sak River and Mae Klong River-include river cross section graphic, river stage, riverbed bottom and the river width, from the Royal Irrigation Department. Bangkok and its vicinity are covered by the Bangkok clay, resulting in no recharging water from the river into the top layer of the model. However, the recharge from the river is still significant to the recharging area outside of Bangkok. The riverbed bottom level such Table 3 of the rivers is laid down on the second layer of the model (BK aquifer). The pumping rate data is for groundwater users authorized by the DGR. The data is estimated by comparing the water usage rate in each province divided by the number of wells. Approximately 8800 wells are distributed over the area shown in Figure 8. Actually, there are still more users that use just a little water but they are not authorized; water smuggling also occurs but there is no data that could be included in this model development.
Water 2018, 10, x; doi: FOR PEER REVIEW 12 of 24 8800 wells are distributed over the area shown in Figure 8. Actually, there are still more users that use just a little water but they are not authorized; water smuggling also occurs but there is no data that could be included in this model development.

Simulation Scenarios
This study models in the steady state by initial calibration with data from 2009 to 2010, since, during this time, the data is quite complete and has little variation, corresponding to the steady flow hypothesis. To obtain a representative hydraulic conductivity value, the model was verified by using the period 2011-2012. The value of the recharge was adjusted every month, including 24 values into the model. The monthly groundwater level of each month was synthesized. Finally, the groundwater level is obtained and compared with the groundwater level forecast by the ARIMAX method.

Unit Root Test
The unit root test was done to consider whether the data is stationary [I(0); integrated of order 0] or nonstationary [I(d); d > 0; integrated of order d], to avoid the spurious regression or the timely nonstationary average and fluctuation data, by using the ADF (Augmented Dicky-Fuller Test) of the groundwater level data from 14 stations.
The initial time series data test found that the groundwater level data was nonstationary. The results from the ADF on values of each station showed that the values were lower than the critical value at the 0.01 and 0.05 significant level and under the null hypothesis acceptance (H0: θ = 0), which means that this data set had unit root. Then, the time series data was adjusted to be stationary by transforming the time series data with the 1st difference. When the adjusted data was retested with the ADF, the results showed that the ADF test statistic values were all higher than the critical value at the 0.01 and 0.05 significance level. This means that the 1st differential data were suitable for developing the ARIMA and ARIMAX model for forecasting the groundwater level in the study area,

Simulation Scenarios
This study models in the steady state by initial calibration with data from 2009 to 2010, since, during this time, the data is quite complete and has little variation, corresponding to the steady flow hypothesis. To obtain a representative hydraulic conductivity value, the model was verified by using the period 2011-2012. The value of the recharge was adjusted every month, including 24 values into the model. The monthly groundwater level of each month was synthesized. Finally, the groundwater level is obtained and compared with the groundwater level forecast by the ARIMAX method.

Unit Root Test
The unit root test was done to consider whether the data is stationary [I(0); integrated of order 0] or nonstationary [I(d); d > 0; integrated of order d], to avoid the spurious regression or the timely nonstationary average and fluctuation data, by using the ADF (Augmented Dicky-Fuller Test) of the groundwater level data from 14 stations.
The initial time series data test found that the groundwater level data was nonstationary. The results from the ADF on values of each station showed that the values were lower than the critical value at the 0.01 and 0.05 significant level and under the null hypothesis acceptance (H 0 : θ = 0), which means that this data set had unit root. Then, the time series data was adjusted to be stationary by transforming the time series data with the 1st difference. When the adjusted data was retested with the ADF, the results showed that the ADF test statistic values were all higher than the critical value at the 0.01 and 0.05 significance level. This means that the 1st differential data were suitable for developing the ARIMA and ARIMAX model for forecasting the groundwater level in the study area, as shown in Table 4.

Identification
Initially, according to the principal of ARIMA modeling, considering the correlation graph of ACF and PACF of the groundwater level in each station, no station showed any sign of nonstationary seasonal form or data which was related to the unit root test result. The correlation graph of ACF showed the long-term regression, while in the correlation graph of PACF, the value quickly lowered to zero at every tested station. The sample of CT4 station in Figure 9 shows that the PACF lowered to 0 in lag 2.

Identification
Initially, according to the principal of ARIMA modeling, considering the correlation graph of ACF and PACF of the groundwater level in each station, no station showed any sign of nonstationary seasonal form or data which was related to the unit root test result. The correlation graph of ACF showed the long-term regression, while in the correlation graph of PACF, the value quickly lowered to zero at every tested station. The sample of CT4 station in Figure 9 shows that the PACF lowered to 0 in lag 2. After adjusting the time series data to be stationary by transforming the initial time series groundwater level data with the 1st difference, or ARIMA(p,1q) and after considering the correlogram of ACF and PACF, such as the time series, the results of the CT4 station are shown in Figure 10. It was found to be a stationary time series. Therefore, the forecasting model is possible, as shown in Table 5. After adjusting the time series data to be stationary by transforming the initial time series groundwater level data with the 1st difference, or ARIMA(p,1q) and after considering the correlogram of ACF and PACF, such as the time series, the results of the CT4 station are shown in Figure 10. It was found to be a stationary time series. Therefore, the forecasting model is possible, as shown in Table 5.  Table 6 shows the possible models of individual study stations, based on the correlogram of ACF and PACF. The model for station CT4 is probably only one form, ARIMA(0,1,2), while station CT5/2 has probably 2 models: ARIMA (1,1,10) and ARIMA(0,2,4). The model that has the lowest value of BIC is ARIMA (1,1,10), which would be used in forecasting for station CT5/2. The other stations are considered in the same way.   Table 6 shows the possible models of individual study stations, based on the correlogram of ACF and PACF. The model for station CT4 is probably only one form, ARIMA(0,1,2), while station CT5/2 has probably 2 models: ARIMA(1,1,10) and ARIMA(0,2,4). The model that has the lowest value of BIC is ARIMA (1,1,10), which would be used in forecasting for station CT5/2. The other stations are considered in the same way. In Table 6, the constant term coefficient is −0.033. The t-statistic value is close to 0, with the 0.05 significance level. This means that the constant term relies on ∆y t . On the other hand, the coefficient of AR(4) is 0.169. The t-statistic value is 3.079, which is far from 0, with the significance level at 0.05. This means that the change of AR(4) is in the same direction with ∆y t . The coefficient of MA(1) is 0.891 and the t-statistic value is 18.618, which is different from 0, with the significance level at 0.05. This means that the change of MA(1) is in the same direction as ∆y t .

Parameter Estimation
From the possible ARIMA model in Table 6, the proper parameters of the model could be found by using the ordinary least square of each station. Table 7 shows the parameter estimation of CT48/2 station, model ARIMA(4,1,1).

Diagnostics
The results of diagnostic checking for white noise consideration of estimated residual (ε t ) show that the correlogram of residuals of autocorrelation (ACF) show no sign of exponential regression. At the same time, the calculated Box and Ljung (Q-statistic) value is lower that the critical value of Chi-square at the 0.10 significance level (prob. < 0.10), which means that ε t is white noise or has normal distribution. The mean is equal to 0 and variances σ 2 , so it could be said that ε t has no autocorrelation and no heteroscedasticity. Table 7 shows the analytic results, showing that all the time series samples passed the diagnostic checking and were suitable for use in forecasting.

ARIMAX
In the ARIMA model, more of the exogenous parameter would be included in the study, such as the climate/oceanology index, IMI, WNPMI, DMI, MEI, SOI and NINO4. All of these were brought in for the Granger Causality Test to see if the index could influence the groundwater level in the study area. The leading parameter test was at the 95% confidence level, so the leading parameter test result was different from station to station. The analytic results of the CT30/1 station are shown in Table 8.  Table 8 shows that the 1-month lag IMI index, 3-month lag MEI index, 3-month lag SOI index and 1-month lag WNPMI index were the leading parameters of the groundwater level. The DMI and Nino index were not the cause of groundwater level at the CT30/1 station, because the test rejected the hypothesis "H 0 : The test index was not the cause of groundwater level" and accepted the hypothesis "H 0 : The groundwater level was not the cause of the text index". The test results of each station for the ARIMAX model are shown in Table 9. Table 9. ARIMAX model of each station.

Station
Model BIC Ljung-Box Q The model from Table 9 can be used to make forecasts in the short time interval for the 14-station monthly groundwater level data in 2011. This was used for accuracy confirmation checking by using the Relative Root Mean Square Error (RRMSE) index. If the RRMSE value was more than 1, the ARIMAX forecasting had less error than the ARIMA model forecasting. Table 10 shows the comparison of RRMSE value between the ARIMAX model and the ARIMA model. At all stations, the forecasting results of ARIMAX model had obviously fewer errors than the ARIMA model, which means that including more climate indices in the ARIMAX model results in more accurate forecasting. After the suitable ARIMAX model had been developed, it was used for the Ex-ante forecast because the ARIMAX model performed very well in the short time interval forecasting. In this study, 12 time intervals were selected for forecasting, January 2012 to December 2012. The forecast value then was compared with the real data. Figure 11 shows the sample of 12-month groundwater level forecasts for CT17/2 station. The forecast value was 19.75 in January 2012 and 19.87 in December 2012. The Root Mean Square Error (RMSE) was 0.719 and MAPE was 3.434% for this station. The results were used for Ex-ante forecasting to compare with the groundwater level data synthesized from the MODFLOW model in the same time interval as the step.  After the suitable ARIMAX model had been developed, it was used for the Ex-ante forecast because the ARIMAX model performed very well in the short time interval forecasting. In this study, 12 time intervals were selected for forecasting, January 2012 to December 2012. The forecast value then was compared with the real data. Figure 11 shows the sample of 12-month groundwater level forecasts for CT17/2 station. The forecast value was 19.75 in January 2012 and 19.87 in December 2012. The Root Mean Square Error (RMSE) was 0.719 and MAPE was 3.434% for this station. The results were used for Ex-ante forecasting to compare with the groundwater level data synthesized from the MODFLOW model in the same time interval as the step.

MODFLOW Simulation in Steady State
After the adjustment of the model using the 325 observed wells distributed in the study area, it was found that the calculated water pressure had the Absolute Residual Mean equal to 7.417 m; the Root Mean Squared Error (RMSE) was 9.489 m; the Normalized RMS was 9.705%; and the Correlation Coefficient was 0.802, as shown in Figure 12.

MODFLOW Simulation in Steady State
After the adjustment of the model using the 325 observed wells distributed in the study area, it was found that the calculated water pressure had the Absolute Residual Mean equal to 7.417 m; the Root Mean Squared Error (RMSE) was 9.489 m; the Normalized RMS was 9.705%; and the Correlation Coefficient was 0.802, as shown in Figure 12. The water level from the model in the steady state has the error in the comparison of the calculated value and the real value measured from the observed wells because there is a large number of observed wells and the study area is vast. Moreover, the assigned hydraulic conductivity for each aquifer had to be the least possible, for the purpose of simplicity of the model. Each aquifer had been set with approximately 2 values-the calibrated vertical hydraulic conductivity and the calibrated horizontal hydraulic conductivity-as shown in Table 11.  The water level from the model in the steady state has the error in the comparison of the calculated value and the real value measured from the observed wells because there is a large number of observed wells and the study area is vast. Moreover, the assigned hydraulic conductivity for each aquifer had to be the least possible, for the purpose of simplicity of the model. Each aquifer had been set with approximately 2 values-the calibrated vertical hydraulic conductivity and the calibrated horizontal hydraulic conductivity-as shown in Table 11. Groundwater flow characteristics of the model assume flow from the edge of the area to the middle of the basin and from north to south. From the retrieved equipotential line and flow direction, the model related to the conceptual model, as shown in Figure 13. From the verification of groundwater level at 325 wells, using the parameter adjustment of the model from the data in 2011, the water level at the time was quite steady, compared with the data from other years. The calculated water pressure value had the Absolute Residual Mean at 9.024 m; Root Mean Squared Error (RMSE) at 11.575 m; Normalized RMS at 11.54%; and Correlation Coefficient at 0.724. The error was seen to be a little higher than those of the calibration.
The sensitivity analysis of the model was carried out by repeating the modeling, by changing the parameters in the model one by one, after the completion of the model adjustment (which related to the conceptual model and had small error), to see the head pressure error caused by changes in each parameter. If any parameter caused a major change, this could mean that that parameter was pertinent to the change in the model.
Changing parameters one at a time and 25% at a time resulted in sensitivity analysis of the model and a better understanding of the influences of the parameters on the model. This information could be used to plan the additional field data collection, to adjust the model for differences in future hydrogeological conditions, or to make the fine parameter adjustments to have smaller errors and to confirm the parameter set correctness before continuing to the next step of model development. From the sensitivity analysis in Figure 14, looking at the relationship between the percentage of parameter changing and the absolute residual mean of head pressure, it was found that the most influential parameter changes were the hydraulic conductivity (K), the pumping rate and RCH value, in that order.
After the process of sensitivity analysis and changing parameters, the verification of the model was repeated with the monthly recharge value at approximately 10% of the monthly rainfall from 2012, in order to synthesize the monthly groundwater level in comparison with the groundwater level retrieved from the ARIMAX model. From the verification of groundwater level at 325 wells, using the parameter adjustment of the model from the data in 2011, the water level at the time was quite steady, compared with the data from other years. The calculated water pressure value had the Absolute Residual Mean at 9.024 m; Root Mean Squared Error (RMSE) at 11.575 m; Normalized RMS at 11.54%; and Correlation Coefficient at 0.724. The error was seen to be a little higher than those of the calibration.
The sensitivity analysis of the model was carried out by repeating the modeling, by changing the parameters in the model one by one, after the completion of the model adjustment (which related to the conceptual model and had small error), to see the head pressure error caused by changes in each parameter. If any parameter caused a major change, this could mean that that parameter was pertinent to the change in the model.
Changing parameters one at a time and 25% at a time resulted in sensitivity analysis of the model and a better understanding of the influences of the parameters on the model. This information could be used to plan the additional field data collection, to adjust the model for differences in future hydrogeological conditions, or to make the fine parameter adjustments to have smaller errors and to confirm the parameter set correctness before continuing to the next step of model development.
From the sensitivity analysis in Figure 14, looking at the relationship between the percentage of parameter changing and the absolute residual mean of head pressure, it was found that the most influential parameter changes were the hydraulic conductivity (K), the pumping rate and RCH value, in that order.
After the process of sensitivity analysis and changing parameters, the verification of the model was repeated with the monthly recharge value at approximately 10% of the monthly rainfall from 2012, in order to synthesize the monthly groundwater level in comparison with the groundwater level retrieved from the ARIMAX model.

Comparison of ARIMAX Groundwater Level and MODFLOW Groundwater Level
The forecast groundwater level from both ARIMAX and MODFLOW were compared by use of the Relative Root Mean Square Error (RRMSE). The results show that the ARIMAX model gave more accurate results than MODFLOW gave. Figure 15 shows the sample of the forecasting at CT17/2 and CT22/3 stations that were done by the MODFLOW model at the 95% confidence interval, which gave results close to the ARIMAX model. In some wells that were outside of the 95% confidence interval, more different forecast results would be presented. It could go up to 100% error in some wells. However, the MODFLOW model analysis in this case was set to be steady flow, so the data brought into the model needed to use the average value, which could have some error. The advantage of using the ARIMAX analysis in this case is the ability to analyze the groundwater level for each well accurately by using only the past groundwater level data and climate indices that were suitable for the case. For MODFLOW, more physical data were needed, such as geographical data, aquifer and river water level, to develop the model. However, there were a number of advantages to using MODFLOW, such as the convenience in importing data about the complex aquifers, ease of data pre-processing and post-processing and the ability to display the regional results, which could clearly show the flow behavior. The selection depended on the analysis objective of each case.

Conclusions
Groundwater is affected both directly and indirectly by climate change. Because groundwater recharge varies with the filling process, groundwater is dependent on many factors, such as the size and type of recharge area and the water recharge rate. These effects can be observed by monitoring

Comparison of ARIMAX Groundwater Level and MODFLOW Groundwater Level
The forecast groundwater level from both ARIMAX and MODFLOW were compared by use of the Relative Root Mean Square Error (RRMSE). The results show that the ARIMAX model gave more accurate results than MODFLOW gave. Figure 15 shows the sample of the forecasting at CT17/2 and CT22/3 stations that were done by the MODFLOW model at the 95% confidence interval, which gave results close to the ARIMAX model. In some wells that were outside of the 95% confidence interval, more different forecast results would be presented. It could go up to 100% error in some wells. However, the MODFLOW model analysis in this case was set to be steady flow, so the data brought into the model needed to use the average value, which could have some error.

Comparison of ARIMAX Groundwater Level and MODFLOW Groundwater Level
The forecast groundwater level from both ARIMAX and MODFLOW were compared by use of the Relative Root Mean Square Error (RRMSE). The results show that the ARIMAX model gave more accurate results than MODFLOW gave. Figure 15 shows the sample of the forecasting at CT17/2 and CT22/3 stations that were done by the MODFLOW model at the 95% confidence interval, which gave results close to the ARIMAX model. In some wells that were outside of the 95% confidence interval, more different forecast results would be presented. It could go up to 100% error in some wells. However, the MODFLOW model analysis in this case was set to be steady flow, so the data brought into the model needed to use the average value, which could have some error. The advantage of using the ARIMAX analysis in this case is the ability to analyze the groundwater level for each well accurately by using only the past groundwater level data and climate indices that were suitable for the case. For MODFLOW, more physical data were needed, such as geographical data, aquifer and river water level, to develop the model. However, there were a number of advantages to using MODFLOW, such as the convenience in importing data about the complex aquifers, ease of data pre-processing and post-processing and the ability to display the regional results, which could clearly show the flow behavior. The selection depended on the analysis objective of each case.

Conclusions
Groundwater is affected both directly and indirectly by climate change. Because groundwater recharge varies with the filling process, groundwater is dependent on many factors, such as the size and type of recharge area and the water recharge rate. These effects can be observed by monitoring The advantage of using the ARIMAX analysis in this case is the ability to analyze the groundwater level for each well accurately by using only the past groundwater level data and climate indices that were suitable for the case. For MODFLOW, more physical data were needed, such as geographical data, aquifer and river water level, to develop the model. However, there were a number of advantages to using MODFLOW, such as the convenience in importing data about the complex aquifers, ease of data pre-processing and post-processing and the ability to display the regional results, which could clearly show the flow behavior. The selection depended on the analysis objective of each case.

Conclusions
Groundwater is affected both directly and indirectly by climate change. Because groundwater recharge varies with the filling process, groundwater is dependent on many factors, such as the size and type of recharge area and the water recharge rate. These effects can be observed by monitoring groundwater levels continuously over a long period. Long-term observations will enable us to plan and manage groundwater use effectively.
The study examined the relationship between the climate index and the groundwater level in the lower Chao Phraya basin by using the climate indices, DMI, IMI, MEI, NINO4, SOI and WNPMI and groundwater level data for 14 stations between 1980 and 2012. In the ARIMAX analysis, the first step was to analyze the stationary characteristics of data before finding the proper model by the Unit Root Test, using the Augmented Dicky-Fuller Test (ADF Test). The results show that the groundwater level was not steady; therefore, the time-series data transformation using the 1st difference was applied, in order to have more appropriate data. The ARIMA model consideration was done by considering the correlogram of ACF and PACF, which gave around 25 possible forms. The BIC statistics were then used to select 1 form per station. The next step was to perform the diagnostic checking to consider the white noise characteristics of the estimated residual (ε t ) by the Box and Ljung process (Q-statistic); in this step, it was found that every selected form was appropriate. The Granger Causality Test of the leading parameters and climate indices was applied to see which index could be used in the ARIMAX model and could forecast the groundwater level for 2012. The results showed that the groundwater level was related to the climate index and could give an effective forecast result at the approximate average RMSE of 0.6. The MODFLOW model then was developed for the study area, with 8 groundwater layers and topped with Bangkok clay. Other boundaries were set to be steady. The verification was done according to the procedure giving the monthly groundwater level result in 2012 at the 95% confidence interval, close to those from the ARIMAX modelling.
The time series forecast of groundwater level by ARIMAX found that it can predict better than ARIMA. It shows that the groundwater level is linked to the climate index and the use of exogenous variables is appropriate. The groundwater simulation model using Visual MODFLOW was found to have advantages in many applications. Importing complicated aquifer information is convenient and it is easy to perform pre-processing and post-processing. However, the definitions of the conditions in this model may not correspond to the actual state of nature, so the simulation results may differ slightly from observation data.
In this study, the MODFLOW model was used with the basic parameter boundary setting and the steady state analysis. The study could be further done in the transient condition. Although, the sensitivity analysis showed that the groundwater pumping was the important parameter, there are still a lot of other unpresented groundwater usage data. This lack of data is one of the cause resulted in the error of the model. If the groundwater pumping influence could be cut out of the groundwater level data [58], the model could then have more chance to get higher accuracy.
Climate change is expected to affect the groundwater resources and the world's populations in the future, so it is going to be a genuine risk factor and will probably become increasingly severe. The study has shown that applying statistical analysis to the study and assessment of the relationship between the climate and the groundwater could be very useful for planning water use for sustainable development. A relatively small amount of data is required and the expense is not very high, either. This could be another way to cope with climate change effectively. In the past few decades, a lot more modeling techniques were presented to the research in hydrological modeling, such as, Artificial Neural Networks (ANNs) [59][60][61]. It would be more useful if these techniques were applied to the groundwater analysis.