Electrical resistivity tomography determines the spatial distribution of clay layer thickness and aquifer vulnerability, Kandal Province, Cambodia

Despite being rich in water resources, many areas of South East Asia face di ﬃ culties in securing clean water supply. This is particularly problematic in regions with a rapidly growing population. In this study, the spatial variability of the thickness of a clay layer, controlling surface – groundwater interactions that a ﬀ ect aquifer vulnerability, was investigated using electrical resistivity tomography (ERT). Data were acquired along two transects, showing signi ﬁ cant di ﬀ erences in the imaged resistivities. Borehole samples were analyzed regarding particle density and composition, and linked to their resistivity. The obtained relationships were used to translate the ﬁ eld electrical resistivities into lithologies. Those revealed considerable variations in the thickness of the clay layer, ranging from 0 m up to 25 m. Geochemical data, highlighting zones of increased ingress of surface water into the groundwater, con ﬁ rmed areas of discontinuities in the clay layer, which act as preferential ﬂ ow paths. The results may guide urban planning of the Phnom Penh city expansion, in order to supply the growing population with safe water. The presented approach of using geophysics to estimate groundwater availability, accessibility, and vulnerability is not only applicable to Kandal Province, Cambodia, but also to many other areas of fast urbanization in South East Asia and beyond.


Introduction
South East Asian countries are naturally rich in water resources. During the monsoon season surface water is typically abundant, while during the dry season water needs to be abstracted from groundwater sources (Babel and Wahid, 2009;Ha et al., 2015). Growing population ( Fig. 1) and food demand is leading to increasing groundwater abstraction, and will further intensify groundwater stress in future years (Erban and Gorelick, 2016;ICEM, 2013). Climate change, which is predicted to cause an increase in temperature of about 10% with negligible change in precipitation (< 5%) in parts of South East Asia, will also have a significant impact on the water availability. It may cause prolonged droughts and changes in the precipitation patterns (ICEM, 2013;Miyan, 2015). While groundwater is abundant and easily accessible through shallow wells, South East Asia, and the Mekong and Red River deltas in particular, are facing challenges in securing clean groundwater resources (Babel and Wahid, 2009); for example, concentrations of arsenic detrimental to human health are well-known in groundwater in these regions (Berg et al., 2001;Bretzler et al., 2017;Charlet and Polya, 2006;Fendorf et al., 2010;Huang et al., 2016;Lado et al., 2008;Polizzotto et al., 2008;Polya et al., 2005Polya et al., , 2003Polya and Lawson, 2015;Winkel et al., 2008). Within the Tonle Sap catchment of Cambodia, an area of growing urbanization, more than 50% of the population are without access to safe drinking water (Babel and Wahid, 2009). In Southern and South Eastern Asia, water scarcity is projected to be one of the major problems faced by urban growth (McDonald et al., 2011).
Over the last 10 years urbanization in South East Asia has resulted in an increase in urban population three times higher than the equivalent increase in the United States (Fig. 1). Densely populated lowland river basins of circum-Himalayan Asia, such as the Mekong River Basin, are facing rapidly growing population rates, and thus further exploitation of already heavily exploited aquifers (Berg et al., 2001;Schneider et al., 2015). The development of mega cities, with inadequate urban planning, is leading to environmental degradation through increased surface run-off and water contamination (Sheng, 2012;Yuen, 2009). Increasing population requires increasing food production, causing a rise in agricultural and irrigated land (Fig. 1).
Although the arsenic in circum-Himalayan groundwater is considered to be of geogenic origin, there is an ongoing debate about the importance of irrigation-enhanced ingress of surface-derived waters (Harvey et al., 2002;Lawson et al., 2016Lawson et al., , 2013Polya and Charlet, 2009). Therefore, an increase in large-scale groundwater withdrawal, for example to meet the demand for irrigation, may further stress the groundwater resource by potentially causing a rise in arsenic concentrations, and may further limit the availability of clean groundwater resources (Aggarwal, 2003;Harvey, 2003;Harvey et al., 2002;McArthur et al., 2011;van Geen, 2003). Kandal Province, Cambodia, is one of several areas in southern and South East Asia that is broadly typical of such arsenic-impacted aquifers (Charlet and Polya, 2006) and is the focus of significant research studies of arsenic biogeochemistry (Charlet and Polya, 2006;Héry et al., 2008;Kocar et al., 2014;Lawson et al., 2013;Polizzotto et al., 2008;Polya, 2008;Polya and Charlet, 2009;Richards et al., 2017;Rowland et al., 2007;Sovann and Polya, 2014). The increasing trend towards urbanization is particularly apparent in northern Kandal Province surrounding Phnom Penh, including a recent expansion of Phnom Penh city limits. The study area forms a major future development zone of the growing Cambodian capital; the population of which is projected to increase by 2030 by more than 75% (ODC, 2008). Thus, it is increasingly important to assess the hydrology and groundwater resources particularly in these rapidly developing regions; to inform future sustainable groundwater management strategies, and to assess groundwater contamination risk in order to secure clean water supply for expanding cities in Cambodia and the greater South East Asian region.
Surfacegroundwater interactions are likely to impact the biogeochemical controls on arsenic release (Richards et al., 2017). While intrusive techniques are able to directly sample the subsurface, they are usually restricted to the sampling of discrete locations, which may be insufficient to fully understand heterogeneous systems and can be logistically very challenging. Geophysical imaging can overcome this limitation by providing spatial or volumetric distributions of subsurface physical properties (Binley et al., 2015;Loke et al., 2013). Geoelectrical techniques, which image the electrical resistivity of the subsurface, have been shown to be able to provide reliable, spatial estimates of surface watergroundwater interactions and groundwater flow characteristics, which are in agreement with conventional intrusive investigations (Cardenas and Markowski, 2011;Chambers et al., 2015;Uhlemann et al., 2016;Williams et al., 2016). This is facilitated by the fact that the resistivity of earth materials is not only a function of their lithology, where clay content in particular usually causes a significant resistivity decrease, but also a function of formation water content and pore water conductivity (where increasing water content and/or salinity usually reduces resistivity; Archie (1942)). Due to the sensitivity of resistivity to lithology, geoelectrical imaging has been employed to produce ground models that are not only comparable to those derived from conventional ground investigation techniques, but offer more detail due to their spatially distributed information (Chambers et al., , 2014a(Chambers et al., , 2014b, thus providing detailed insights into the heterogeneity of sedimentary deposits (Parsekian et al., 2015;Van Dam, 2012).
Geophysical imaging, and electrical resistivity tomography (ERT) in particular, is inherently non-unique, meaning that the measured data are explained equally well by an infinite number of subsurface physical models, and has a spatially variable sensitivity (Binley et al., 2015;Loke et al., 2013;Olayinka and Yaramanci, 2000). Thus, geophysical imaging is preferably interpreted together with some kind of intrusive investigation (Chambers et al., 2014a(Chambers et al., , 2014bComas et al., 2004;Merritt et al., 2013;Salas-Romero et al., 2015), which can significantly reduce the non-uniqueness and increases the reliability of subsurface models and interpretations.
In order to understand the lithological controls on arsenic groundwater contamination in northern Kandal Province, Cambodia, geoelectrical measurements have been performed to image the spatial extent and continuity of a subsurface clay layer, which is expected to act as a hydraulic constraint on surface watergroundwater interactions. By combining field resistivity measurements with laboratory analysis of particle size and resistivity of borehole samples and field intrusive investigations, we present a detailed picture of the variability of this "protective" layer. This study also provides a substantial addition to the sedimentological characterization of the study area.

Study area
This study focuses on an area of approximately 50 km 2 within northern Kandal Province, Cambodia, located approximately 12 km southeast of Phnom Penh in the lower Mekong River Basin (Fig. 2). South-East of Phnom Penh, the Mekong River divides into the Bassac River and a continuation of the main Mekong River, which form the major distributaries of the Mekong Delta. The latter can be subdivided into two parts: an upper delta plain which is dominated by fluvial processes and within which the study area is located, and a lower delta plain that is characterized by marine processes (Nguyen et al., 2000). The Mekong delta was formed within the last 8 ka, prograding more than 200 km, largely in response to high sediment loads from the Mekong River, and Quaternary sea-level changes (Tamura et al., 2009(Tamura et al., , 2007. While the Holocene sediments of the lower delta plain are well studied (Nguyen et al., 2000;Ta et al., 2002;Xue et al., 2014), detailed geological investigations of the deposits of the upper delta plain are largely limited to those reported by Tamura et al. (2007Tamura et al. ( , 2009. A 30.7 m drill core recovered and analyzed by Tamura et al. (2007), from a borehole located less than 10 km south-east of the present study area, has been interpreted as a sequence of fluvial and tidal depositional facies, and is summarized in ascending order as follows; (a) fluvial, (b) tidal, (c) salt marsh, and (d) flood-plain deposits, with (e) natural-levee to flood-plain deposits forming the upper 3-4 m. The (e) natural-levee to flood-plain deposits show a reddish brown sandy silt and high mud content. The underlying (d) flood-plain deposits comprise a layer of greenish gray clays, around 3 m thick, which shows a sharp lower boundary with the underlying (c) salt marsh deposits. The latter are characterized by a 5 m thick alternation of peat and graded sand to silt. The peat comprises laminated gray clay to silt and shows centimeter sized plant fragments, with sandy layers formed of yellowish sand and silt. The boundary to the underlying (b) tidal deposits is gradational. Around 16 m in thickness, the tidal deposits consist of very fine sand with mud drapes, and a mud content that generally decreases with depth. Together with the cross-laminated micaceous sands of the underlying (a) fluvial deposits, they form the aquifer in this area. This aquifer is confined by the poorly permeable clays of the (c) salt marsh and the (d) flood-plain deposits. In previous studies, this clay-rich layer was found to be continuous, with thicknesses ranging between 5 m and 22 m . Due to different mechanisms of electrical conduction, clay deposits generally show significantly lower resistivities than sand deposits. Thus, for this study the (c) salt marsh and the (d) flood-plain deposits interpreted as an aquitard were expected to show lower resistivities than and a sharp contrast to the underlying (b) tidal and (a) fluvial sand deposits forming the aquifer.
Water in this area is mainly produced through shallow domestic wells (commonly around 30 m in depth), supplying houses and farms with drinking and irrigation water. Usually, these are drilled to a depth only a few meters below the clay layer. The water quality is known to be problematic, having high arsenic and manganese concentrations (Buschmann et al., 2007;Feldman et al., 2007;Huang et al., 2016;Lawson et al., 2013;Polizzotto et al., 2008;Polya et al., 2005;Richards et al., 2017;Rowland et al., 2008;Sovann and Polya, 2014). The main industry in this region is agriculture, with predominately small scale rice and fruit plantations (ICEM, 2013). The presence of infrastructure is very limited (Fig. 2).

Geoelectrical imaging
Geoelectrical data were collected in May 2013, towards the end of the dry season. Data were acquired along 12 ERT profiles, ranging in length between 83 and 230 m. These profiles were located along two transects (Fig. 2), usually next to local roads and tracks. Retrospectively, they were named "T-Sand" and "T-Clay" in accordance with their main respective lithologies (Richards et al., 2017), which were determined by the resistivity measurements. An electrode spacing of 2 m was used for all profiles except for the short profile P24, which was only 83 m long and where a 1 m electrode spacing was used (Table 1). Resistivity measurements are usually conducted using four-point electrode arrangements; two electrodes to introduce current into the ground and two electrodes to measure the resulting potential. The measurement sequence comprised conventional dipole-dipole measurements, using typical dipole lengths a of 2-14 m (1-7 m for profile P24, and 2-12 m for profiles P14 and P23) and dipole separations na with n = 1-8. The dipole-dipole array was chosen as it is known to provide high resolution data while effectively exploiting multichannel measurement systems (Chambers et al., 2002;Dahlin and Zhou, 2004). This choice was supported by numerical simulation of two plausible geological scenarios (Fig. 3). The first case represents a resistive (200 Ω m, e.g. sandy) overlying a more conductive (10 Ω m, e.g. clay) layer. The second case, somewhat more complex, comprises a thin, discontinuous resistive layer with conductive elements, overlying a more conductive layer, which in turn overlies a resistive layer. Both models include a thin (0.5 m) surficial layer of intermediate resistivity (75 Ω m). From those models synthetic data were calculated for standard dipole-dipole (a = 2-12 m, n = 1-8) and Wenner configurations (a = 2-36 m). Those synthetic data were then inverted using the same parameters as described for the field data set in the following paragraph. For the first case, dipole-dipole and Wenner yield comparable results, both being able to image the interface between the conductive and underlying resistive layer at the correct depth. For the second case, the dipole-dipole survey images depth and structure of the shallow resistive and underlying conductive layer to a higher accuracy than the Wenner survey, which overestimates the depth of the conductive layer and introduces some spurious structure. At depth, the Wenner survey images the resistivity of the deeper model parts more accurately than the dipole-dipole survey. LaBrecque et al. (1996), Olayinka and Yaramanci (2000), Tsourlos et al. (1999) and Chambers et al. (2002) present more detailed studies on the effects of noise, model complexity, and array type on the inversion results.
The measurement system employed in this study was the AGI SuperSting R8/IP instrument. Each resistivity measurement was made in normal and reciprocal configuration (LaBrecque et al., 1996), where the reciprocal measurement is equal to the normal measurements but with interchanged injection and measurement dipoles, with the measured value being defined as the mean of these two measurements. The measurement error was defined as the standard error in the mean; this is commonly referred to as the reciprocal error in ERT.
Data quality was generally good to very good, with less than 13.3% of the data having reciprocal errors of more than 5%. A threshold of 5% reciprocal error was also used for data filtering; all measurements exceeding this value were removed from the data set. The average measurement error ranged from 0.18% to 1.16%, with the highest errors being found at P25, which was situated along a track between houses of a village. To obtain a resistivity model of the subsurface, the measured resistances need to be inverted. ERT inversion is fundamentally a nonunique problem in that an infinite number of models can explain the measured data equally well (Olayinka and Yaramanci, 2000). Thus, to define an optimum solution the inversion process needs to be constrained. Here, the data were inverted using a smoothness-constrained least-squares inversion method, employing a robust data weighting and a L2-norm for the model roughness (Loke and Barker, 1996). The sensitivities of the Jacobian matrix were calculated using a full Gauss-Newton approach. Most inversion routines use a "global" damping factor. However, recent work has shown that using a spatially variable damping factor defined by the sensitivity may provide superior results (Yi et al., 2003). This approach, known as "active constraint balancing", was employed in this study. The geological setting was expected to comprise mostly (sub-)horizontal layers, thus the inversion was constrained by a horizontal flatness filter, placing twice the weight on horizontally continuous features than vertical. The electrode positions and local topography were determined using dGPS equipment (C-Nav 3050 Dual Frequency & Glonass GPS receiver, C-Nav, USA, with Star-Fire wide-area differential GPS corrections) and included in the inversion process by employing a finite-element mesh. In the inversion, data were weighted by their measured reciprocal error, provided it was larger than a lower limit of 10% of the average noise level. This was calculated from the error values between the 10% and 90% percentiles . A commonly used measure of goodness for the inversion is the root-mean-squared (RMS) misfit between modelled and measured resistances. For all profiles, the inversion converged in less than 7 iterations to small RMS misfits between 0.4 and 1.7% (Table 1), highlighting the good agreement between measured data and inverted resistivity models. For data visualization, the Depth-of-Investigation Index (DOI) was employed (Oldenburg and Li, 1999). Regions of the resistivity models with DOI's of more than 0.2 were displayed semitransparent, as for these regions the resistivity structures are less well determined by the measured data.

Intrusive investigation and sampling
Boreholes were installed in transects perpendicular to the Mekong and/or Bassac Rivers using local manual rotary drilling techniques as described by Richards et al. (2015). Borehole locations were targeted to be as close as possible to ERT profiles (see Fig. 2); however, in some cases this was practically restricted due to space availability and/or local landowner permissions. Borehole locations can be divided into two types, major and minor, where major locations consisted of five to six nested boreholes targeting depths between 6 m and 45 m, while minor locations consisted of two nested boreholes at 15 and 30 m depth. Wet sediment cores were collected at the time of drilling using a locally-designed stainless steel sampler, fitted internally with a replaceable extruded acrylic tube (25 or 50 mm outer diameter) capped on the bottom with a replaceable one-way stopper for sample collection. The sampler was screwed onto 3 m segments of narrow drilling rod and inserted into the drilling pipe until the desired sampling depth was reached. Sediment samples were collected by manual hammering of the sampler into the ground until sufficient sample was collected and retrieved from the sampler immediately upon collection. Sediment samples were targeted every 3 m of drilling depth, as practically feasible, up to a maximum of 45 m. Sediment subsamples were stored separately for (i) particle size analysis and resistivity measurements (frozen in polyethylene bags) and (ii) X-ray fluorescence spectrometry (XRF; frozen under nitrogen in pre-furnaced aluminium foil packets stored in polyethylene bags).

Laboratory analysis
Particle size distribution was measured at the British Geological Survey (Keyworth, UK) on 110 dried and sieved (< 2 mm) sediment samples using a laser diffraction particle size analyzer (LS13 320, Beckman Coulter) enabled with polarization intensity differential scattering as previously described (Richards et al., 2017). Samples were sonicated for five minutes to break down clumped particles prior to particle size determination (run time 90 s). A certified control (Coulter LS G35D Garnet, Beckman Coulter) was analyzed approximately every 20 samples to verify the instrument accuracy. A subset of seven samples was run in triplicate and precision was estimated to be generally within 1% for each size fraction.
The particle-size fractions are defined as (i) clay < 8 µm; (ii) silt 8 µm -0.063 mm; and (iii) sand 0.063 -2 mm, reflecting a deviation from the conventional clay/silt size cutoff of 2 µm; this was adapted because of the challenges sometimes associated with measuring the size of non-spherical particles with laser-based methods (Rawlins et al., 2009). Sedimentary classification has been done using an adapted Shepard's ternary diagram based on the relative proportions of sand, silt and clay sized particles (Shepard, 1954). Given the deltaic (fluvial/ marine) depositional environment, the now widely accepted textural definition of mud being all sediments finer than 0.063 mm (e.g. silt and clay combined) has also been used (Flemming, 2000).
Samples designated for elemental analysis using X-ray fluorescence spectrometry (XRF) were prepared by drying and grinding samples into a fine powder using a grinding mill (Planetary Mono Mill Pulverisette 6, Fritsch). Ground sample (12 g) was mixed with a crystalline wax (3 g, Hoechst Wax C Micropowder) and pressed into pellets using a hydraulic press. Major and trace elemental composition of pelleted sediment samples was measured using XRF (Axios Sequential XRF Spectrometer, PANalytical) at the Manchester Analytical Geochemistry Unit (University of Manchester, UK). Twelve sediment samples were pelleted and analyzed in triplicate and errors in XRF were estimated to be within 2%.
Laboratory resistivity measurements on a subset of 17 sediment samples were conducted at the British Geological Survey for comparison to modelled resistivity outputs from ERT field measurements. The subset covered the range of particle sizes and sample depths. Wet sediments were compacted into small, sealed boxes of dimensions 25 mm × 25 mm × 75 mm. Current electrodes were placed at the narrow ends of the boxes using stainless steel plates extending across the full width and height of the box, ensuring a homogeneous current flow through the sample (Hen-Jones et al., 2017). Potential electrodes were placed in the centre of the sample with a spacing of 20 mm. Current was injected using a current source, generating a square wave of 10 Hz with an amplitude of 2.7 mA. The potential drop across the two potential electrodes was measured using a digital oscilloscope (Tektronix TDS2024B). The sample resistivity is then given as: with the measured voltage drop V, the transmitter current I tx , and the geometric factor K. The geometric factor can be calculated analytically by: with the sample cross-section A and the spacing between potential electrodes l. This equation assumes a homogenous current flow through the sample. The actual design of the sample holder and dimensions of the sample were thought to influence the effective value of K; therefore the sample holder and different sample sizes were modelled in Comsol Multiphysics® (COMSOL, 2011) and the results verified using liquids of known resistivity. Through drying and wetting of the samples, resistivity was determined at moisture contents of 15, 20, 25, and 30%, which are representative of the field conditions.

Resistivity models and borehole logs -T-Sand
T-Sand consisted of six ERT profiles, sampling a total distance of 1060 m up to 30 m depth. These profiles were located along a field track between the Bassac and Mekong Rivers in a broad SW to NE direction. Distances to the Bassac ranged from 150 to 4550 m, and to the Mekong between 2450 and 6750 m. Resistivities imaged along this transect range between 4 and 460 Ω m, with a median of 28 Ω m, and show two to three distinct horizontal layers. Generally, lower resistivities were imaged close to the Bassac than further inland. Profile P11, located only about 150 m NE of the Bassac River, shows a thin resistive top layer (Fig. 4e), reflecting the compacted field track, overlying low resistivities with fine layers of higher resistivity, which are dipping toward the Bassac, throughout the section. Further inland (Profiles P14 and P15, Fig. 4d and c), resistivities close to the surface are similar to those defining the majority of the profile close to the Bassac. Within this upper layer horizontal features of increased, but also decreased resistivity can be found, indicating fine-grained deposits with higher sand or clay content, respectively. This overlies a significantly more resistive layer. Profile P16, located furthest away from the Bassac is missing the low-resistivity upper layer and shows high resistivities throughout the profile, interrupted by small-scale low-resistivity features, some of which are showing horizontal alignment (Fig. 4b).
To investigate the subsurface in a transverse orientation to this broad SW-NE direction, P12 and P13 were acquired roughly perpendicular to each other (Fig. 4f), with P12 being oriented in a NW-SE direction. Both profiles show a similar vertical sequence with a ∼10 m thick layer of the lowest resistivities imaged along this transect, sandwiched between two more resistive layers. Within this low-resistivity layer, horizontal features of elevated resistivity can be found, forming an intermittent thin layer at less than 5 m below ground level (bgl). While most structures are spatially consistent between these two adjacent profiles, two resistive features at about 15 m bgl are evident in P13 only, and may represent spatially confined sand lenses. Fig. 4 also shows the logs of boreholes located close to the ERT profiles. For T-Sand, most of the boreholes show a similar, spatially consistent sequence, with a shallow horizon of clayey silt overlying a > 25 m thick layer of sand to silty sand. This layer is underlain by another clayey facies, which is fining with depth from a clayey silt to a silty clay formation. This deep clayey facies can be associated to Pleistocene deposits that are known to be of high mud content (Tamura et al., 2009(Tamura et al., , 2007. Thus, the borehole logs correlate well with the imaged resistivities, which show a similar stratification.

Resistivity models and borehole logs -T-Clay
The resistivities imaged along T-Clay range between 3.5 and 1080 Ω m, with a median of 17 Ω m, and are generally lower than those of T-Sand. Six profiles were acquired along this transect, located 500-2750 m south of the Mekong river, sampling a distance of 913 m up to a depth of about 25 m bgl. Profiles P22 and P23 were acquired along field tracks, P21 and P24 are situated in open grassland, while P25 and P26 where acquired next to and between houses and roads. Profile P26 (Fig. 5b), located closest to the Mekong River, shows a transition from low-resistivity to highly resistive material with increasing distance to the Mekong. This highly resistive material forms a wedge, which thins out towards the Mekong and shows a thickness of more than 15 m at the south-western end. Profiles P25 and P24 show a two-layer sequence of a less than 5 m thick, resistive surface layer, overlying a significantly more conductive layer, which extends below the depth of investigation of these profiles. At P25 (Fig. 5c), this layer shows small-scale, spatially delimited features of elevated resistivity, which are likely to be regions dominated by coarser sedimentary grain size. This is in agreement with a borehole log acquired at this site, showing a sequence of silty clay, sandwiched between layers of clayey silt. These small-scale features are missing at P24, which shows a homogeneous, conductive lower layer, with increasing resistivities towards greater depths. These findings correlate well with samples taken from a nearby borehole.
A similar sequence can be observed at profiles P21 to P23, which were acquired close to each other with P21 and P22 oriented in a NW-SE direction, and P23 in a SW-NE direction (Fig. 5e). Here, the lowresistivity layer merely has a thickness of less than 10 m, and overlies a more resistive layer. The boundary between these two layers is gradational. This observation is in agreement with a borehole located between P22 and P23, which shows gradational coarsening from a clayey facies over silty clay and clayey silt to a silty sand and sandy facies, at depths comparable to those of lithological boundaries derived from the resistivity images.

Laboratory ground-truthing
When plotted on an adapted Shepard's diagram (Fig. 6a, Shepard (1954)), the particle size analysis classified 51% of the obtained samples as silty sand, 17% as clayey silt, and 15% as sand. The remaining 17% classified as silty clay (9%), clayey sand (5%), sandy silt (2%), and sandy clay (1%). Using Folk's nomenclature (Folk, 1968), nearly all samples fall within the categories of sand, muddy sand, sandy mud, and mud, highlighting their high contents of fines and their sedimentary environment within a flood plain. Samples taken from shallow depths above 10 m bgl showed generally higher mud content (clay and silt) than samples taken from depths below 20 m bgl, which had higher sand content (Fig. 6a). This is in agreement with the local geology and depositional environment (Tamura et al., 2007), where (d) flood-plain and (c) salt-marsh deposits with very high to high mud content are overlying (b) sandy tidal deposits of significantly lower mud content. In a nearby borehole, the boundary between them was found at a depth of 12 m bgl, which is comparable to our results. Fig. 6a also indicates the smallest distance of the samples to one of the two rivers. The sediment classes are distributed over all distances. Although a large number of samples classified as sand originate less than 2000 m from the river, and a high proportion classified as silty sand can be found at more than 4000 m from the rivers, the distance is nevertheless unlikely to be a major control on sediment class.
Resistivity of the borehole samples was analyzed with respect to variations in sediment composition and moisture content (Fig. 6b). The data show the expected decreasing resistivity with increasing clay and silt content (Archie, 1942;Waxman and Smits, 1968), which is caused by a change in the contribution of the two electric conduction principles in sediments, electrolytic and surface conduction. While the resistivity of clean sand is mainly determined by the electrolytic conduction through the pore fluid, clay minerals also allow for a surface conduction through an electric double layer (Telford et al., 1990). Thus, increasing clay content causes a decrease in bulk resistivity. The data highlight that sediment texture (as reflected in the particle size distribution) is the controlling factor of the sample resistivity, while variations in moisture content have a limited impact. For all three main sediment classes, sand, silt, and clay, decreasing moisture content leads to an increase in resistivity, which is less significant than the resistivity variations caused by changes in grain size. Thus, the resistivity data were fitted against sand, silt, and clay content to allow for a translation of field resistivity measurements to sediment classes. Silt and clay fractions were fitted to the sample resistivity using exponential functions, and sand content using a linear function (Fig. 6b). While the data is well explained by the fit for sand and clay content (Pearson's r = 0.61 and r = 0.67, respectively), for silt a slightly greater spread can be observed (Pearson's r = 0.59). This is likely to be caused by the transition between the main electrical conduction mechanisms, which has a stronger control on the resistivity of sand and clay, respectively. These functions were used to determine sand, silt, and clay fractions from the field resistivities.
Comparative analysis was carried out between elemental compositions of the samples for silica, aluminium(III) and iron(III) oxide obtained from the XRF analysis and sample resistivity, sand, and clay fraction (Fig. 7). This analysis showed a positive linear correlation between silica fraction and sample resistivity, as well as between silica and sand fraction. As silica is an electrical insulator and the main constituent of sand, these correlations were expected. Increasing aluminium(III) oxide content relates to increasing clay content and decreasing resistivity. Elevated aluminium oxide content is typical for highly weathered tropical soils rich in clay minerals (Schaetzl and Anderson, 2005). No correlation between iron(III) oxide content and resistivity or sand and clay fraction was found.   6. (a) Classification of sediments according to particle size distribution determined from samples taken at depths between 3 m and 45 m, and at varying distances from the river (modified after Shepard (1954)). (b) Resistivity measurements plotted against sand, silt, and clay fraction of a subset of 17 samples from the data presented in (a). Resistivity measurements were taken at four different moisture contents, 15%, 20%, 25%, and 30%. Fig. 7. Cross-plots of XRF results with resistivity, sand, clay fraction; samples of T-Sand are shown as dots, and samples of T-Clay as diamonds. With increasing silica content, resistivity increases, while increasing aluminium(III) oxide relates to a decreasing resistivity. These trends are independent of the sample moisture content. Silica and sand content, and aluminium (III) oxide and clay content show a positive linear correlation as expected. No correlation between iron(III) oxide and resistivity was found.

Resistivity-derived particle size imaging
Particle sizes derived from the ERT measurements were translated into sedimentary classes following Shepard's nomenclature using thresholds described in Smith et al. (2009). The results are shown in Figs. 8 and 9 for T-Sand and T-Clay, respectively. The distribution obtained from all ERT profiles was comparable to the laboratory particle size analysis, with sand and silty sand accounting for about 20% each, clayey silt for about 25%, silty clay 15%, sandy silt for about 10%, and silt and clay for less than 5% each. Note that this translation has some intrinsic limitations, namely the non-uniqueness and spatially variable sensitivity of ERT (Day-Lewis, 2005), where, particularly with increasing depth, the resistivity and thus the derived lithology are becoming less reliable.
Similarly to the resistivity results of T-Sand, profiles close to the Bassac River (P11 to P13) show high clay and silt content, with profiles P12 and P13 having a ∼10 m thick clayey layer, which is overlain by a less than 2 m thick surficial layer of silty sand to sand (Fig. 8). Beneath the clayey facies another sandy layer is imaged, which can be interpreted as (b) tidal deposits forming the local aquifer. Further inland, the uppermost sandy layer of natural levee deposits is missing and the clayrich (d) flood-plain and (c) salt marsh deposits, being up to 7 m thick, are reaching the surface. However, it is shown to be an intermittent feature (P14 and P15), forming clay-free "windows" of potentially high permeability, thus allowing for surfacegroundwater interaction. At P16, located centrally between the Bassac and Mekong River, the clayey facies is either less than 2 m thick or entirely absent. The sandy tidal deposits are thus exposed at the surface, hence it is likely that surface waters can easily infiltrate into the groundwater. The resistivity-derived lithology is in good agreement with the borehole data, which are consistently showing a thick sandy layer that is overlain in places by a thin surficial layer of clayey silt.
T-Clay (Fig. 9) shows significantly higher clay and silt contents. All profiles, with the exception of P26, show a thick clayey silt to clay facies, which is overlain by a surficial, thin layer of silty sand, and underlain by a sandy layer, forming the local aquifer. This agrees with the sequence observed by Tamura et al. (2007) of sandy (e) natural-levee to flood-plain deposits, overlying clay-rich (d) flood-plain and (c) salt marsh deposits, which in turn overlie the sandy (b) tidal deposits. The clayey facies ranges in thickness between 15 m and possibly more than 25 m. Its high clay content is likely to significantly reduce its permeability, thereby providing a continuous barrier for surface water infiltration into the underlying aquifer. Within profiles P21 to P24, it is shown to be reasonably homogeneous, with increasing grain size at depth, and almost horizontally aligned clay layers. At P25, 1500 m from the Mekong River, this clay layer becomes patchy and intermittent, with less distinct horizontal layering, and significantly more heterogeneity, interrupted by isolated lenses with increased grain size. Profile P26, located about 500 m from the Mekong, is characterized by a thick sandy layer that thins towards the Mekong, and contains isolated clay lenses. A borehole next to this line confirmed the high sand content, showing silty sand and sand down to 30 m bgl. The other boreholes

Discussion
The resistivity-derived lithologies and borehole logs were used to investigate the spatial variability of the thickness of the low-permeability clay layer overlying the sandy aquifer. The thickness of this shallow clay layer (comprising units classified as clay, silty clay, and clayey silt) was extracted from the ERT profiles at 5 m intervals and, together with the borehole data, used as input for an interpolation of the shallow clay thickness over the wider study area. Data were interpolated onto a grid with 25 × 25 m cell size using kriging, by employing an exponential variogram model with a range of 450 m, a sill of 30 m, and a nugget effect of 2 m. The results show a clear difference between T-Sand and T-Clay. Along T-Sand, the thickness of the surficial clay layer varies between 0 and 15 m, while at T-Clay profiles P21 to P25 are located in an area with clay thickness of more than 15 m. Within the otherwise thick clay layer of T-Clay, profile P26 is located in a spatially localised area where the clayey facies are either thin or entirely absent. This shows that the clayey flood-plain and salt marsh deposits, which protect the underlying sandy aquifer from rapid surface water ingress, vary significantly in thickness between completely absent to more than 25 m, thus creating potential "hotspots" of increased surface watergroundwater interaction. This is in contrast to previous studies in this region Polizzotto et al., 2008), which suggested a continuous, 5-20 m thick surficial clay layer protecting the aquifer from ingress of surface water.
This sedimentological interpretation is consistent with observations from aqueous geochemical data obtained from the same study area. In particular, the stable isotopes of hydrogen and oxygen have been used to provenance groundwater recharge sources and possible interactions between various water bodies in this system using mixing models (Richards et al., 2017). These models suggested that recharge in the zones identified in Figs. 8-10 as having very thin surface clay coverage also exhibited a stable isotopic signature similar to local precipitation, consistent with the hypothesis that these are zones of rapid recharge to the highly permeable (b) tidal sands. Conversely, the groundwater in areas shown to have thicker clay layers in Fig. 10 exhibits more evaporative isotopic signatures, which may be indicative of different recharge processes in these locations. Furthermore, the inorganic geochemical characterization of groundwater shows seasonal changes in groundwater chemistry (including key parameters like sulfate and dissolved oxygen), which are likely associated with monsoonal patterns driving rapid surface-groundwater interactions in certain zones of the study area (Richards et al., 2017).
The accuracy of the clay thickness map depends mainly upon two factors, the accuracy of the resistivity model and the validity of the translation from resistivity to lithology. As described already, resistivity imaging is inherently non-unique and sensitivity decreases with increasing depth. Thus, deep parts of the resistivity model are less reliable than shallow areas. In this study we tried to address this limitation by presenting the data that corresponds to areas of small sensitivities (DOI > 0.2) transparent, as those regions of the models are more unreliable. Another limitation of resistivity imaging are the smoothness constraints applied in the inversion, which, even if sharp resistivity contrasts exist in the subsurface, will cause the inversion to produce smooth subsurface resistivity models. Thus, the derived lithological models will also present a "smoothed" image of the true subsurface lithology.
Resistivity is largely controlled by mineralogy, whereby an increase in the fraction of clay minerals usually lowers the resistivity, but pore water saturation, conductivity, and material temperature also have an effect. Temperature is likely to decrease with depth and would cause a smooth increase in resistivity (Chambers et al., 2014a(Chambers et al., , 2014bHayley et al., 2007). In comparison to other factors influencing the subsurface resistivity, this effect is deemed negligible. Water saturation levels and pore water conductivity are known to influence bulk subsurface resistivities (Archie, 1942;Telford et al., 1990), whereby decreasing saturation can cause a significant increase in resistivity. In this study, the effect of varying saturation on the subsurface resistivity was tested in laboratory experiments (Fig. 6b) for a range of moisture contents representative for the field conditions during the survey; sample resistivities were determined at 15-30% gravimetric moisture content. The most significant resistivity variation was found for samples with a high sand fraction (> 60% sand), where a change from 30% to 15% moisture content caused a rise in resistivity of up to 100%. For samples with a low sand fraction, and thus high silt and clay fraction, this variation was significantly smaller and mostly below 30%. The high variability for samples with a high sand fraction was expected, as for sandy materials the main mechanism for current conduction is electrolytic. Comparing samples with small (< 20%) and high (> 80%) sand content shows a resistivity range between 7.5 Ω m and 39 Ω m at the same moisture content (30%), which clearly outweighs the variation caused by changes in moisture content.
While the relationship between resistivity and sand and clay content is relatively close, the spread is wider for the resistivity-silt relationship. This may introduce another potential source of error for the translation of resistivity to lithology. As we are mainly interested in the "end" members, i.e., sand and clay, this higher uncertainty for the silt-derived lithologies is acceptable and should not have a significant effect on the estimation of the surficial clay thickness, as our definition of the clay layer comprised all clayey lithologies. The good agreement between lithologies derived from the resistivity measurements and logged from the boreholes, supports the reliability of the resistivity models and the validity of our approach. It shows that combining geophysical with intrusive and laboratory studies can significantly improve the reliability of the derived subsurface models.
This study, despite focusing on only two transects, shows that the thickness of the clay layer, controlling surfacegroundwater interactions, is highly variable in the study area. Linking this observation with the geochemical investigations revealed that the most vulnerable areas of the underlying sandy aquifer are those where the clay layer is missing or very thin, thus providing potential pathways where organic rich surface waters could rapidly infiltrate. Despite those areas offering a seemingly attractive opportunity to produce water from shallow wells, the water quality is likely to be problematic due to potentially increased arsenic and other contaminant concentrations. In order to reach less vulnerable parts of the aquifer, locations are more favourable where a sufficiently continuous, poorly permeable clay layer protects the aquifer. These conclusions may inform urban planning in the context of the Phnom Penh city development, in terms of where best to locate water production facilities that will inevitably be needed. However, in our attempts to develop a more comprehensive understanding of the variability in clay thickness and to assess the impact of permeable windows acting as preferential flow paths for organics-rich surface waters into the groundwater, this study can only form an initial step. To cover a wider area in sufficient detail, additional resistivity transects and airborne surveys based on electromagnetic methods could be used; these can also measure the subsurface resistivity and are known to be able to determine lithologies even in complex environments (Parsekian et al., 2015;Pryet et al., 2011;Robinson et al., 2010). The aquifer vulnerability, but also its production potential, is a function of the hydraulic conductivity of the different units, which can be estimated through pumping or tracer tests, but also inferred from monitoring of resistivity variations (Chambers et al., 2015;Kemna et al., Fig. 10. Map of clay thickness derived from the resistivity derived lithology and borehole logs. The data were interpolated using kriging, by employing an exponential variogram model with a range of 450 m, a sill of 30 m, and a nugget effect of 2 m. 2002; Kuras et al., 2016;Slater et al., 2002;Uhlemann et al., 2017). Focusing such efforts on those permeable windows allowing for a fast infiltration of surface waters into the groundwater may provide insight into flow process characteristics, and when combined with geochemical analysis may further enhance our understanding of arsenic release mechanisms.
The approach presented here is not only applicable to the groundwater problems and geological setting of the wider Phnom Penh area, but also to most expanding cities that are in need of safe water supply. Being non-intrusive and fast, geophysical investigations can provide data at a spatial density and coverage that is not achievable using conventional intrusive testing and sampling. Thus, they can provide supporting data on the local lithology that may inform groundwater availability, accessibility, and vulnerability; factors that are crucial for city development and for addressing water shortage (McDonald et al., 2011).

Conclusions
Access to safe drinking water is becoming increasingly difficult in areas of fast urbanization, particularly in developing countries. In this study, we employed geoelectrical measurements to investigate the Holocene deltaic deposits controlling the hydrology of Kandal Province, Cambodia, an important area for the expansion of Cambodia's capital, Phnom Penh. Geoelectrical data were acquired along two transects with the aim of assessing the variability of a clay layer that is thought to control surfacegroundwater interactions, thereby affecting the aquifer's vulnerability. Resistivity data showed significant differences between these two transects, with "T-Sand" exhibiting considerably higher resistivities than "T-Clay". Through linking laboratory analysis of particle size distributions and XRF measurements with sample resistivity, the field resistivities were translated into lithologies. While transect "T-Clay" showed an up to 25 m thick, mostly continuous layer of clay protecting the aquifer from ingress of organic rich surface waters, at "T-Sand" this layer was either absent or less than 15 m thick. The clay layer thickness was extracted from the resistivity profiles and, together with borehole information, interpolated onto a wider area. It was shown to be highly variable, ranging from 0 m up to 25 m. Complementary geochemical data not only confirmed locations of permeable windows through this clay layer, but also highlighted the increased aquifer vulnerability where the protective clay layer is absent or discontinuous.
Understanding the local hydrogeological settings, controlling water availability, accessibility and vulnerability, is vital for planning of urban developments in order to provide the expanding population with safe drinking water. For this, the presented approach is not only applicable to Kandal Province, Cambodia, but also to many other regions within South East Asia and beyond, which face rapidly growing urban populations and potential threats to clean water supply.