3D stratigraphic modelling of the Bangkok basin using Kriging on borehole data

The Lower Central Plain of Thailand has a deep and highly irregular basement filled with complex layers of alluvial sandy soil and deltaic clay or silt. The Bangkok Metropolitan region with its high level of infrastructure development is located in this plain. With high population density, the problem of land subsidence is critical. This study uses borehole data and the Kriging method to interpolate the data. A detailed 3D stratigraphic model of the basin is presented and several cross sections along two directions. Seventeen layers and the points used for modelling each layer are presented. The model shows eight aquifers in the Bangkok basin, lying between eight layers of clay, and a bottom layer above the basement. The bottom of the lowest aquifer of the Bangkok aquifer system is 610 m deep. The basement extends down to a depth of more than 2000 m in some places. Our 3D model, which extends to the basement, is of interest for mining and site-specific seismic risk analysis. Moreover, our results can be very useful for groundwater and land subsidence studies.


Introduction
The central region of the Kingdom of Thailand has seen intense economic development and construction growth. The areas in and around Bangkok are the regional hub of transport and logistics, with high level of development [1]. The Bangkok Metropolitan Region covers an area of 7,762 km 2 . The planning region includes Krung Thep Mahanakhon, commonly known as Bangkok and five adjacent provinces: Nakhon Pathom, Nonthaburi, Pathum Thani, Samut Prakan and Samut Sakhon. The region has a growing populationmore than 10 million people in 2000, more than 14 million in 2010, and earlier projections estimated a population of 16 million for 2019. The population density in the region is high -5631 people per km 2 in Krung Thep Mahanakhon [2]. Excessive deep well pumping has led to critical land subsidence in Bangkok [3]. To carry out aquifer studies and implement ground-water resource management, a 3D stratified model and knowledge of the geotechnical properties of the local strata is necessary. For regions with high infrastructure development, such as Bangkok, such information is important in many fields and applications. The mechanical properties of the soil or rock is paramount in geotechnical engineering when selecting the type, depth and size of structural foundations suitable for the bearing capacity of the soil. It is also necessary for landform system studies, for establishing the groundwater table and for soil or water contamination studies. The stratigraphy and geotechnical data are also significant for mining, construction of subsurface structures and material surveys. For seismic risk estimation and risk reduction, it is necessary to study site-specific amplification of seismic waves, which requires knowledge of the local stratigraphic model and the geotechnical properties of the area.
In 2014, 3D geological modelling and a study of the geotechnical characteristics of the subsoils was conducted in Phnom Penh City, Cambodia [4]. In a similar study, a 3D geological model was developed for Bangkok subsoils, to a depth of up to 60 m [5]. Most of the subsurface surveys carried out in the Bangkok basin were for land subsidence studies. Although shallow subsurface profiles of Bangkok have been established, there is lack of information on deeper soil stratification and a detailed layer profile of the basin has not been published. This study uses the available information for deeper strata to construct a stratigraphic model for the deeper layers of Bangkok subsoils.

Geological Setting
The Lower Central Plain or Chao Phraya plain ( Figure 1) is a large flat low land plain [3]. The approximate width of the Chao Phraya plain in the North-South direction from Chainat to the estuary of the Chao-Phraya river in the Gulf of Thailand is 200 km and the total area of the plain is approximately 36,000 km 2 . The plain is mildly sloped with an average elevation of 2 m above sea level [6].
Data from boreholes show that the Chao Phraya basement is large, with a depth of more than 2000 m, and composed of granites and metamorphic rocks. The basement depth varies considerably in the different regions of the basin, including near the estuary of the Chao-Phraya river. In some areas, the depth of basement is unknown. The basement is filled with interbedded alluvial and flood-plain silt of the Chao-Phraya delta. The sequence of aquifers and clays is known only down to 450 m depth. The sequence has eight confined aquifers. The top layer of the sequence is soft marine clay called Bangkok clay, covering an area of 13,800 km 2 of the plain [3].

Data and Methodology
Borehole data from a previous study on groundwater flow [7] was used in our study. The data from 113 boreholes ( Figure 1) show 17 layers of alternating clay strata and sand aquifers. The 17 layers are listed in descending order in Table 1. While not all the boreholes reached the bottom layer, many were between 100 and 200 m depth. The borehole points for each layer are shown in Figure 2. The properties of layers BC to PN of the Bangkok basin are well established in other literatures. The borehole data used in our study has one more aquifer layer 'T' below layer PN.
To construct the 3D-model, the data for each layer was interpolated between the gridded points in the study area. Adding depths in each layer for every grid point, a 17-layered profile was created. The result estimated from the Spline method passes through the input values, but estimates may exceed the range of the input data. Kriging is a geostatistical method, which uses correlation between the data rather than a fixed set of mathematical equations and is an exact interpolator. For dense and uniformly distributed data locations, all the methods provide good and similar estimated results. For clustered data locations with large gaps, the values cannot be well estimated by any of the methods. However, Kriging gives less weight to the clustered points or treats clusters as single points. This helps to compensate for the effect of clustering and gives better estimates than the other interpolation techniques. Kriging also calculates the estimation error at the estimation points [8]. There are various Kriging methods based on the trend of the data. Simple Kriging (SK) assumes a constant trend with a known mean. Ordinary Kriging (OK) assumes a constant mean in the neighbourhood of the estimation point, rather than the whole domain. Universal Kriging (UK) or Kriging with a trend is a variant of OK where a linear or higher order trend is fitted instead of a local mean in the neighbourhood of the estimation point [8]. In Kriging, the empirical semivariogram calculated is replaced with an acceptable semivariogram model. Choosing the semivariogram model is based on subjective judgement. Spherical and exponential models are appropriate for representing properties with higher levels of short-term variability [8]. For a lag distance 'h', range 'a' and sill 'c', a spherical model can be represented as Equation (1) and an exponential model can be represented by Equation (2): The borehole data in this study is not uniformly distributed within the study area. There is a large cluster near the estuary of Chao-Phraya basin in the Gulf of Bangkok while in surrounding locations the boreholes are sparsely distributed (Figure 1). The surface of the basin usually does not have good correlation between values separated by large distances. In other words, the variable may have a spatial trend. To analyse such a dataset, the UK method is the most suitable of all the methods mentioned above. We also analysed the data using the IDW, NN and Spline methods. For the clustered dataset, IDW did not give satisfactory estimates. The estimates of the Spline method were far higher than the input values. The estimates of the NN and UK methods were better than those of the IDW and Spline methods. The surfaces produced by the UK method were smoother than those of the NN method. Based on the theoretical suitability as well as visual fitting, UK was used to interpolate all seventeen layers.
In the UK method, a trend surface is fitted and residuals from the trend are calculated. We used the 'gstat' package in R [9]. An empirical variogram is generated from the data and fitted into an existing acceptable semivariogram model. The gstat package has a choice of exponential, spherical, Gaussian or Matern models. In this study, either the exponential or spherical model were used based on the nature of the empirical variogram and the least-weighted sum of the squared errors of the fitted model (SSErr). The range coefficients and partial sill coefficients were fitted manually or automatically by the program, using the default fitting method.    Figure 2 shows the interpolated values of each layer. The red points are the borehole data points for each layer. The blue points represent the points with no sediment overlying the bedrock, based on a previous study [10], which were added to the dataset to confine the region of interpolation. The points are dense in the upper layers of the aquifer system; the point density gradually decreases towards the lower layers, especially in the centre of the study area at the Ang Thong and Phra Nakhon Si Ayutthaya provinces. For the bottom layer, there are no points available in those provinces ( Figure 2). However, in the Bangkok metropolitan region data is available for all the layers. Layers NB and NL extend through a large area of the basin whereas layer T is confined to a small area.

Results and Discussion
The 2D-cross sections along six lines in the East-West direction and six lines along the North-South direction are shown in Figure 3(a) and (b), respectively, from layer BC to PN of the Bangkok aquifer system. The maximum depth of the bottom of the last aquifer layer of the Bangkok aquifer system is 610 m. Research on site-specific seismic amplification requires information of the deeper layers. Figure  3(c) shows the cross-sections including layer T for three lines along the E-W and four lines along the N-S directions. The lines are shown in Figure 1. The depth of the aquifer system decreases towards the north of the basin.
The maximum basement depth is 2109 m at [13.546N, 100.477E]. At one point in the Singburi province, the depth is 1550 m. Sixteen data locations are available for layer T, which is not dense enough in the centre of the study area. The cross section shows a considerable decrease in depth for layer T from 14º to 14.75º in section Y4 of Figure 3(c). The results cannot be verified for layer T in the Singburi region; however, in the area surrounding Krung Thep Mahanakhon, where more data is available, the results are reliable for all the layers.

Conclusion
In this study we built a 3D stratigraphic model of the Bangkok basin using borehole data. The model shows 17 layers of varying thickness. The individual layers and the points used for modelling the layers are presented. The Bangkok aquifer system comprises eight aquifers confined by alternating clay layers. The maximum depth of the bottom of the lowest aquifer of the Bangkok aquifer system is 610 m. The Bangkok basin is more than 2000 m deep in places. Cross sections along two directions were produced for clear representation of the results.
The accuracy of such analysis depends on the density of the data. The data is denser in the area surrounding the Krung Thep Mahanakhon, where the model can be used with confidence. The 3D stratigraphic model presented here can have applications in broad fields such as groundwater studies, land subsidence studies, and site-specific seismic risk analysis.