Characterization of karst structures using quasi-3D electrical resistivity tomography

Karst is characteristically complex, hydrogeologically, due to a high degree of heterogeneity, which is often typified by specific features, for example, cavities and sinkholes, embedded in a landscape with significant spatial variability of weathering. Characterization of such heterogeneity is challenging with conventional hydrogeological methods, however, geophysical tools offer the potential to gain insight into key features that control the hydrological function of a karst aquifer. Electrical resistivity tomography (ERT) is recognized as the most effective technique for mapping karstic features. This method is typically carried out along transects to reveal 2D models of resistivity variability. However, karstic systems are rarely 2D in nature. In this study, ERT is employed in valley and hillslope regions of a karst critical zone observatory (Chenqi watershed, Guizhou province, China), using a quasi-3D approach. The results from the extensive geophysical surveys show that there is a strong association between resistivity anomalies and known karstic features. They highlight the significance of a marlstone layer in channeling spring flow in the catchment and confining deeper groundwater flow, evidenced by, for example, localized artesian conditions in observation wells. Our results highlight the need to analyze and interpret geophysical data in a three-dimensional manner in such highly heterogeneous karstic environments, and the value of combining geological and hydrogeological data with geophysical models to help improve our understanding of the hydrological function of a karst system.


Introduction
Karst accounts for 7-12% of the world's land area (Ford and Williams 2007). Karstic aquifers supply almost a quarter of the world's population (Ford and Williams 2007;Hartmann et al. 2014;Stevanovic 2018) and, in many regions are vulnerable to land use and, consequently, changes in land management. Karst aquifers are notoriously complex, including three types of porosities, i.e., micropores, small fissures and fractures, and large fractures and conduits, which result in significant multi-scale heterogeneity (Bakalowicz 2005). Water storage dynamics and the movement of water and solutes in karst systems are subject to duality because of these physical features of karst (Hartmann et al. 2014). To develop conceptual and, subsequently, numerical models of the hydrology of a karst aquifer it is necessary to assess the nature and significance of such features.
Assessing aquifer features in karstic environments using traditional drilling approaches is limited because of the high degree of heterogeneity (Goldscheider and Drew 2007). Similarly, interpretation of conventional pumping tests and borehole water hydrochemistry is challenging. Geophysical methods can provide some insight into the heterogeneity of karst aquifers, and ultimately be used in combination with conventional borehole-based methods to develop conceptual models of aquifer function. Electrical resistivity tomography (ERT) (e.g., Gunther et al. 2006;Binley 2015) is widely used in groundwater investigations and is recognized as an effective technique for revealing karst aquifer structure (e.g., Al-Fares et al. 2002;Chalikakis et al. 2011;Margiotta et al. 2012Margiotta et al. , 2016Metwaly and Al Fouzan 2013;Kaufmann and Deceuster 2014). Furthermore, time-lapse ERT imaging has been successfully used to assess groundwater pathways in karstic environments (e.g., Sawyer et al. 2015;Watlet et al. 2018).
Electrical resistivity tomography provides 2D or 3D images of the variation of electrical resistivity of the subsurface using measurements made with electrodes, typically on the ground surface and/or in boreholes. The electrical resistivity of geological materials is a function of lithological properties (e.g., pore size distribution), degree of saturation, and pore water composition (e.g., Lesmes and Friedman 2005). Generally, in karstic aquifers, a water-filled void (e.g., fracture, conduit and cavity) has a lower resistivity than the surrounding rocks (e.g., limestone and dolomite), whereas an air-filled void has a higher resistivity than the host rock (Zhu et al. 2011). The interpretation of electrical resistivity can, however, be challenging. For example, the existence of clay filled voids or clay rich bedrock can lead to low resistivity features due to both high pore water retention and elevated electrical conductivity of grain boundaries. Therefore, low resistivity features are not necessarily associated with high hydraulic conductivity. The magnitude of resistivity contrasts can be site-specific due to other contributing factors (e.g., Robert et al. 2011;Xu et al. 2016;Bermejo et al. 2017;McCormack et al. 2017;Keshavarzi et al. 2017), highlighting a limitation of using electrical geophysics in isolation of other observations. Electrical resistivity tomography is typically used to derive 2D images of the subsurface, i.e., it is assumed that resistivity does not vary in the strike direction (normal to the ERT survey line). Whilst this is often plausible for horizontally layered systems (e.g., Cardarelli and De Donno 2017), in karstic aquifers, resistivity values are likely to be highly variable in three directions. Therefore, 2D models may not correctly reflect the resistivity anomalies, both in terms of size and position (Chávez et al. 2015). Kaufmann and Deceuster (2014) highlighted that 3D ERT surveys should be preferred to 2D ERT, especially in areas where the main directions of fractures and similar features are unknown.
Most previous electrical geophysical studies have focused on using ERT to detect single features of a karst system, such as sinkholes, cavities or conduits (e.g., Chalikakis et al. 2011;Meyerhoff et al. 2012;Gutierrez et al. 2014;Parise et al. 2018). Studies of entire watersheds are rare, and Hartmann et al. (2014) points out that the availability of data limits the characterization of karst aquifer heterogeneity. In this study we report on an extensive geophysical field campaign in a small research watershed of the karst critical zone observatory in Guizhou province, China. Our focus is the improved characterization of the hydrogeological structure of the karst environment using electrical geophysics, with a particular focus on the use of quasi-3D geophysical modeling techniques.

Site description
The study site is the Chenqi watershed located in Puding County, Guizhou Province, China (Fig. 1). The Chenqi watershed is a sub-catchment of the Houzhai watershed. It has an area of 0.92 km 2 , and a classical cockpit karst landform ( Fig. 1, 2). Chenqi has been the focus of numerous studies of karst hydrology (e.g., Zhang et al. 2013). The average elevation in the catchment is approximately 1387 m a.s.l. with the standard deviation of 46.5 m, and mean slope of 24°. The lower slopes of most hillslopes in the catchment have been transformed into terraced fields by local farmers (see Fig. 2), which possibly benefits from the layered rock formations with a gentle dip angle. In the main valley (as illustrated in Fig. 2) most of the area has been cultivated for paddy fields. The climate belongs to the subtropical humid monsoonal region. Average annual precipitation and temperature are 1315 mm and 15.1 °C, respectively.
The geological sections ( Fig. 1) demonstrate that there are four kinds of exposed rocks: thick limestone, dolomite, lamellar limestone (i.e., thin-bedded limestone) and marlstone. The thick limestone and some dolomite are located in the upper layers, i.e., hilltops, which are the youngest rocks in the Chenqi watershed. Outcrops of lamellar limestone occur on hillsides, and, in contrast, the interbedded limestone and marlstone can be seen in the bottom of hillslopes. Figure 3 shows an example exposed profile of limestone and marlstone in the catchment (location of the photograph is shown in Fig. 1).
In the Chenqi catchment, there are no obvious faults, and the rock exhibits a near horizontal layered sedimentary structure (Figs. 1, 3). The average inclination direction and dip angle in the geological sections ( Fig. 1) are approximately 270° and 7°, respectively. The rock layers are generally inclined towards the catchment outlet, which affect the flow direction and the karst landform evolution. The geological sections in Fig. 1 show evidence of earlier connectivity of exposed units (see, for example, the clear connectivity between the main lithological units in section B-B′ and C-C′ in adjacent hillslopes), which have subsequently become disconnected through dissolution. Rock cores extracted within the watershed valley indicate a combination of unfractured limestone and fractured limestones, with some clay infill of fractures.
The main soil type of the Chenqi watershed is a brown clay formed by carbonate after its dissolution in a hot, humid and rainy climate condition (Zhou et al. 2012). Because of serious soil losses (from erosion) in the catchment, the soil layer in the hillslopes is typically shallow The photograph was taken from the east of the watershed; view is westerly and discontinuous (Peng and Wang 2012;Hu et al. 2015;Fu et al. 2016). The soil layer in the valley is much thicker (typically greater than 1 m thick). Therefore, the hillslopes show exposed karst rock, whereas in the valleys the karst rock underlies substantial soil cover. Beneath the surface or soil exposed at the surface, there is an epikarst layer that consist of highly weathered soluble rock (Williams 2008). In our study area, the thickness of the epikarst layer is estimated to be approximately 10 m and the average porosity is about 5% (Zhang et al. 2013).
In the main valley a surface channel exists (blue line in Fig. 1). Water flows in the channel during the wet season. During the dry season the channel is dry, but is used to channel pumped groundwater for irrigation downstream. There are two naturally formed sinkholes in the main valley of the watershed (④ and ⑤ in Fig. 1), which are categorized as bedrock collapse sinkholes (Waltham et al. 2005;Parise and Gunn 2007;Gutierrez et al. 2014). These sinkholes are connected to the main surface channel. An outlet spring (① in Fig. 1) and two hillslope springs (② and ③ in Fig. 1) are natural outlets of groundwater. For a long term study of the groundwater hydraulic dynamics, five wells were previously drilled in the watershed valley (⑥-⑩ in Fig. 1). It is interesting to note that well ⑦ is an artesian well, of which the discharge is steady, about 83 m 3 /day; well ⑧ (15 m from well ⑦), in contrast, is not artesian.
In a recent study, Chen et al. (2018) utilized hydrometric data alongside measurements of water isotopes in rain water and groundwater to partition event and pre-event water. Their study, was supplemented with a small number of 2D electrical resistivity profiles. Their analysis of tracer data illustrates the complex dynamics of groundwater flow in the catchment and support a conceptual model of fast and slow groundwater flow zones in the area.

ERT
Electrical resistivity tomography uses measurements of the potential field, created by injecting current into the ground, to determine the spatial variability of electrical resistivity. It is, therefore, an 'active method'. Current is injected using a pair of electrodes and the potential field is sampled using combinations of pairs of other electrodes. By injecting current in multiple current dipole configurations it is possible to 'probe' different regions of investigation. Increasing the spacing between the current electrodes and/or the distance between the current and the potential dipole, allows a greater depth of investigation (Binley 2015). In this study we limit the electrode spread to target a depth of investigation of up to 30 m.
Given a set of multiple four electrode measurements that sense different areas of the subsurface, inverse methods are applied to determine the 'best' resistivity model that is consistent with the data. As the problem is non-linear, least squares methods based on Gauss-Newton search techniques are commonly employed enhanced with spatial regularization to constrain the inversion process. The iterative process runs until a solution is obtained that matches the observed data to a satisfactory level of misfit. For more information on the details of the process see, for example, Binley (2015).
Typically, such four electrode measurements are made along a transect on the ground surface. Assuming no variation in resistivity normal to the direction of the transect; it is possible to determine a 2D resistivity model that is consistent with the observed data, i.e., a 2D inverse model is developed. However, as stated earlier, in highly heterogeneous environments such assumptions of two-dimensionality may be invalid. In such cases 3D modeling is necessary. Ideally ERT measurements would be made in 3D arrays (i.e., using a plane of electrodes, rather than a single line), however, such a configuration is often impractical because: (1) access to parts of a site may be restricted; (2) cable and electrode requirements for 3D surveys limit surveys to small areas; (3) ERT systems are limited to addressing up to (typically) 100 electrodes at a time. However, it is possible to conduct multiple 2D surveys and combine these in 3D analysis, i.e., develop a quasi 3D model. In this study, we use the 3D inversion program R3t (Binley 2013) which is based on an Occam's solution to the 3D inverse problem, and employs an unstructured tetrahedral finite element mesh for forward modeling and resistivity model parameterization. R3t has been widely used in other (nonkarst) studies (e.g., Koestel et al. 2008;Perri et al. 2012).
Electrical resistivity tomography measurements can be made with different quadrupole geometries. The Wenner, Schlumberger, and dipole-dipole configurations are the most popular (e.g., Kaufmann and Deceuster 2014;Binley 2015;Keshavarzi et al. 2017). Although the dipole-dipole configuration has, in comparison to others, a weak signal strength, it is the most effective for assessing lateral variation in resistivity (a likely characteristic of karst). Zhou et al. (2002) pointed out that the dipole-dipole array is the most effective configuration for mapping sinkholes. In addition, many modern ERT measurement devices allow some level of multi-channel measurement using a dipole-dipole configuration, thus making this configuration efficient in the field. For the surveys reported here we used Syscal Pro 96 (Iris Instruments, France), which is capable of measuring up to ten channels (dipoles) simultaneously.
The locations of ERT survey lines are shown in Fig. 1. The survey areas were selected to target specific features within the catchment to help understand more about the subsurface structure and hydrological functioning. All surveys were carried out during April-June 2017. The ERT lines are located in four areas: main valley channel; (A) hillslope, which includes two hillslope springs (② and ③) and two wells (⑨ and ⑩); (B) middle of the valley, which includes the artesian well ⑦, observation well ⑧ and sinkhole ⑤; (C) outlet, which includes the outlet spring ①, observation well ⑥ and sinkhole ④.
A single ~ 550 m long ERT survey was carried out along the main valley channel as an initial reconnaissance to assess the extent of resistivity variability in the region covering springs and sinkholes, and thus help design subsequent targeted surveys. The survey used 5 m spaced electrodes. The maximum spacing between current and potential dipoles was 50 m, given a depth of investigation of ~ 30 m. The survey included a full set of reciprocal measurements (e.g., Tso et al. 2017) for data error analysis. Data quality was good (reciprocal errors typically less than 2%).
The hillslope spring group was targeted to help understand the hydrological function of the two springs (② and ③) and contrasting water level responses in wells (⑨ and ⑩). 5 ERT lines were surveyed to cover the area (Fig. 1). A total of 8900 ERT measurements were collected, again with a full set of reciprocals for error analysis. A 3D unstructured finite element mesh consisting of approximately 653,000 elements was created to combine the 5 ERT lines and honor the local topography. The local 3D finite element mesh and the location of all 309 electrodes are shown in Fig. 4.
A total of 12 ERT lines were surveyed to target the artesian well group (B) (see Fig. 1 for location of the 12 lines). The group includes a sinkhole ⑤, an artesian well ⑦ (20 m deep) and a 22.3 m deep observation well ⑧ (Fig. 1). The zone was targeted to assess why wells ⑦ and ⑧ show completely different behavior despite their close proximity to each other. The combined ERT dataset consists of 6644 measurements and a full set of reciprocal data for error analysis. Again, a 3D unstructured finite element mesh was created (with approximately 713,000 elements).
Group C contains a sinkhole ④, an observation well ⑥ and an outlet spring ① (Fig. 1). The zone was selected to assess whether there is evidence of connectivity between the spring and sinkhole in this area. 4 ERT lines (Fig. 1) with a total of 6379 measurements (plus a full set of reciprocal measurements) were combined. The local refined grid approach was used again for mesh generation, resulting in a mesh with approximately 447,000 elements.
For the main valley channel we applied the 2D ERT modeling code (Binley 2016). For the other three areas we utilized R3t (Binley 2013). For 3D modeling and inversion, mesh generation in irregular geometry (including topography) can be challenging. Furthermore, fine mesh discretization is required close to electrodes (where potential gradients are high) and thus for large arrays inversions may be limited by available computer resources. R3t allows the use of an unstructured tetrahedral finite elements mesh to help alleviate such challenges. The mesh generator Gmsh (Geuzaine and Remacle 2009) was used to construct all 3D finite element meshes in this study.

Resistivity along the main valley channel
The resistivity model from main valley channel ERT survey is shown in Fig. 5. The resistivity varies between ~ 20 Ωm and ~ 6000 Ωm. Generally, in karstic aquifers, low resistivities can be associated with the high clay content (e.g., marlstone), or soils and weathered with a relatively high water content (Zhu et al. 2011). Following previous studies (Robert et al. 2011;Xu et al. 2016;Bermejo et al. 2017;Keshavarzi et al. 2017) and an understanding of the characteristics of Chenqi watershed (e.g., brown clay and layered rock stratum), the resistivity distribution can be divided into three categories: (i) a soil layer, including Quaternary deposits (resistivity < 100 Ωm); (ii) extensively weathered rock (100 Ωm < resistivity < 400 Ωm); (iii) compact limestone (resistivity > 400 Ωm). Although this may be somewhat simplistic, it helps identify common features (at least based on their geoelectrical properties). Figure 5 shows that the upper low resistivity zone becomes more pronounced in an up-valley direction, suggesting a thickening of soil and unconsolidated deposits. The high resistive regions are localized, as expected in such a karstic environment. An extensive zone of mid-resistivity is seen at depth in the upper valley region (x > 150 m in Fig. 5), which may indicate more permeable shallow bedrock. This low resistivity zone may be related to the openness of valley and impact of erosion-deposition events from historic major floods.
The results from the 2D ERT survey along the main channel was followed by more targeted quasi 3D investigations of the specific nature of the resistivity structure in regions of the valley, as described below.

Hillslope spring group (A)
Historic water level time series and water isotope sampling has revealed contrasting behavior of wells ⑨ and ⑩. Well ⑨ shows damped water level fluctuation that is consistent with unconfined conditions, whereas well ⑩ reveals more confined behavior from water level time series (Chen et al. 2018;see also, Zhang et al. 2013 for analysis of spring response). Similarly from water isotope analysis (Chen et al. 2018) water in samples from the two springs and well ⑨ appears relatively young, in contrast to those from well ⑩.
Combining the 5 ERT surveys in a 3D inversion results in the model shown in Fig. 6. The image shows sections of the model aligned with the 5 transects, where sensitivity of the data is high. Three main features stand out from the image in Fig. 6. (1) In most areas the resistivity is high suggesting extensive unweathered bedrock, with outcrops close to the ground surface, particularly upslope. This high resistivity zone also extends to the area adjacent to well ⑩ (see insert Fig. 6b), consistent with the confined behavior previously noted.
(2) Close to well ⑨, intermediate resistivities suggest extensive weathering/fracturing at depth. This is consistent with the observed unconfined nature of the water levels in the well. (3) Adjacent to the two springs (and in a similar position in the most northerly transect) a low resistivity zone exists. These low resistivity features are wedge-shaped, which is further illustrated by the isosurface from the 3D model in Fig. 7, which indicates that the spatial extent of the feature is about 50-80 m upslope from the springs. It would appear that these wedge features are connected, although we recognize that there is limited off-transect sensitivity of the ERT data. This suggests that there may not be a single

Artesian well group (B)
Inversion of the data resulted in the image shown in Fig. 8. The image shows an extensive zone of high resistivity at depth, but note that adjacent to the artesian well a low resistivity zone exists at depth. In this case we attribute the low resistivity to extensive fracturing and water filled secondary porosity. It would appear that this very localized feature maps the hydraulically conductive pathway that results in artesian conditions, and does not extend 15 m away to the non-artesian well ⑧. This suggests that a pathway for deep upwelling is extremely localized, although the potential for other upwelling sources may exist (based on the other zones of low resistivity mapped at depth in Fig. 8). The 3D image also shows that the sinkhole ⑤ exists in a low resistivity zone that is laterally extensive, with a typical depth of about 8 m, perhaps indicating a shallow unconfined region above the lower permeability limestone. Figure 9 shows the resulting resistivity image from 3D inversion of the combined dataset. It can be seen that the low resistivity area around the outlet spring ① and sinkhole ④ is shallow. In contrast, the low resistivity area close to the observation well ⑥ is much thicker. In comparison to the region further up the valley the unweathered bedrock appears much closer to the ground surface. Note also in Fig. 9 the low resistivity wedge feature at the bottom of the hillslope in the most easterly line-this has similar features to the spring zone noted in group (A).

Location of marlstone layer
The inverted resistivity models in Figs. 6 and 9 reveal low resistivity wedge features near the foot of two hillslopes. From field observations, a marlstone layer was noted along the road at the foot of the southerly hill (Figs. 1, 3). Given that marlstone is typically 35-65% clay and 65-35% carbonate (Pettijohn 1975), we expect this to have a relatively low resistivity (e.g., Martínez-Moreno et al. 2015; Al  Dulaymi et al. 2012). Furthermore, any hydraulic impedance of water flow from above will reduce the resistivity. Therefore, we can deduce that the wedge feature is a marlstone layer. Given the geology profiles and know dip angle of 3° and inclination direction of 270° of strata (Fig. 1), along with a digital terrain map, it is possible to create a 3D representation of the catchment geology. Figure 10 shows such a representation along with the outcrop of a marlstone layer. The road at the foot of the hill (where the photograph in Fig. 2 was taken) is also shown in Fig. 10, which independently aligns with the modeled outcrop.
Examination of Fig. 10 reveals that the marlstone layer in the hillslope spring zone (A) is likely to be the same with that around the hill-foot road. Similarly, the low resistivity wedge feature in the group (C) image from the outlet (Fig. 9) is also likely to be caused by the same geological feature. It appears that the hillslope springs exist as a consequence of the extensive marlstone layering-water originating from Fig. 8 Inverted resistivity distribution for well group (marked (B) in Fig. 1)   Fig. 9 Inverted resistivity distribution for outlet group (marked (C) in Fig. 1) recharge upslope is impeded vertically and, under relatively wet conditions, discharges are seen in the springs. This supports the previous observations of relatively young water emerging in the springs (Chen et al. 2018).
To evaluate the effect of marlstone layer on ERT signals, a resistivity model for the hillslope spring area was created, using the same mesh as before (Fig. 4). In this model a 5-m thick 50 Ωm laterally extensive layer was inserted in a 2500 Ωm background. Using the same measurement scheme for the 5 ERT lines collected in the field, a forward response was calculated with R3t. These data were then inverted to assess what shape feature would result. Figure 11 shows the resulting model. In this figure, the inverted resistivity distribution is similar with that in Fig. 6, adding further evidence that the low resistivity wedge features are, in fact, likely to be a result of a laterally extensive hydraulically impeding unit. The exposed area of the spring zone in Figs. 6, 7 is larger than that in Fig. 11, possibly because of the assumed thickness of 5 m being somewhat greater than the true conditions. Driven by the field ERT survey results and supported by the synthetic modeling results, we can conclude that there is a marlstone layer in the foot of the hillslopes of the Chenqi watershed, and that this layer may be the cause of spring flow in the catchment.

Flow features of springs and sinkholes
In the Chenqi watershed, 2 sinkholes, 2 hillslope springs and 1 outlet spring (Fig. 1) have been found. These spring and sinkholes are distributed in the vicinity of the central axis of the watershed, and are connected by the channel (Fig. 1), which means that during the high rainfall season, not only springs but also sinkholes can act as discharge points.
The resistivity values around the hillslope springs (Figs. 6, 7) all shows that springs exist at transition points in resistivity. A low resistivity unit, associated with the marlstone and pore water retention, appears to act as a horizontally continuous impeding layer which results in the discharge of groundwater at the surface. Figures 6, 7 reveal the springs at the base of the low resistivity unit, which we attribute to the accumulation of pore water above the thin marlstone layer and possible enhanced weathering of the limestone above it. A schematic diagram of the spring feature is shown in Fig. 12.
Similarly, the resistivity values close to the sinkhole ④ ( Fig. 9) are high, and the soil and weathered rocked layers are thin or non-existent. However, upstream of sinkhole ④, the soil and weathered rocked layers are much thicker (Fig. 9). Therefore, it appears that there is a migration of subsurface water to sinkhole ④. When groundwater flow is great, water emerges out at the surface at some point that becomes a broken point in the flow line. Under such conditions this broken point gradually evolves into a sinkhole (e.g., sinkhole ④).
The flow feature of sinkholes is similar to that of springs in the Chenqi watershed (Fig. 12). The difference is that at a spring point, the downward flow is prevented, which results in continuous discharges at the spring(s). In contrast, at a sinkhole, when groundwater flow is low, water still permeates into the ground, but when groundwater flow is great, the water emerges out at the surface. A similar structure can be also found around the outlet spring ① and sinkhole ⑤ in the resistivity images (Figs. 8,9).

Other characteristics of the karstic aquifer
In the hillslope areas surveyed, the soil layer is relatively thin and discontinuous (Fig. 6), which is in agreement with previous studies (Peng and Wang 2012;Hu et al. 2015;Fu et al. 2016) that have pointed out that the soil thickness in hilltops and hillsides are about 20-40 cm and 50-150 cm, respectively. As the permeability of marlstone layer will be low and the unit extensive, the amount of water recharging to the deep aquifer from the hillslope may be limited. Most of water stored in the hillslope is likely to flow out in the form of spring water, which supports isotopic observations: the age of spring water ② (average δD value of − 65) is newer than that of all water bodies (e.g., wells and sinkholes) in the valley (δD of − 60).
In the valley, the thickness of the shallow karstic aquifer is non-uniform and in general appears to gradually decrease from upstream to downstream (Fig. 5). The resistivity surveys, supported by observations from boreholes, also reveal that a deeper confined system is likely to exist in the area. The artesian well ⑦ appears to connect to such a system and the resistivity surveys reveal the localized nature of variability in this region of the catchment. The result is confirmed by geochemical observations from water samples from the well: the artesian well water contains high concentration of sulfate ions (SO 4 2− ), from gypsum, whereas in other (unconfined) boreholes, the concentration of sulfate ions is much lower.

Conclusions
Karst aquifers are complex systems with high heterogeneity, making conceptualization of flow pathways challenging using only traditional borehole-based methods. In the Chenqi catchment such complexity has been observed for some time using such methods. ERT surveys reported here have helped improve our understanding of the hydrological processes in the catchment. A strong association between resistivity features and known karstic structures (such as springs and sinkholes) and contrasting unconfined and confined behavior of wells, has resulted in greater insight into the significance and these features on a catchment scale. This study did not focus on mapping conduits in the Chenqi watershed. Such a focus in the future could draw support from time-lapse ERT (e.g., Meyerhoff et al. 2012).
Traditional 2D resistivity surveys can be effective in karstic environments but if significant lateral variability exists, one needs to consider better the 3D nature of the system under study, which has been the focus here. 3D modeling tools were used to help overcome the challenges of traditional ERT methods. Our surveys were not fully 3D and such surveys are unlikely to be a practical option at scales above a few tens of meters, however, by combining multiple 2D data in a 3D inversion we can, at least, account for some lateral variability away from the survey lines. And where survey lines intersect, a consistent model of resistivity results (which may not be the case for independent 2D imaging).
The study has revealed the geoelectrical contrasts associated with a marlstone layer, and further investigation has highlighted that this marlstone layer may be extensive and be the primary cause of hillslope springs, and limit deeper aquifer recharge across the catchment. The geophysical surveys help explain previously observed behavior of environmental tracers. The resistivity images have also help develop hypotheses for the development of sinkholes within the catchment, i.e., contrasting lateral permeability leading to ephemeral surface discharge (spring flow) under high rainfall conditions, with subsequent collapse under normal drier conditions. Geophysics should never be used in isolation. Coupled with extensive measurements using traditional borehole-based methods, they can, however, offer immense data value, particularly in such complex subsurface environments.