Quantitative estimations of aquifer properties from resistivity in the Bolivian highlands

Groundwater management and aquifer assessment rely on the reasonable knowledge of its hydrogeological characteristics, i.e., transmissivity (T), hydraulic conductivity (K), storativity (S) and specific storage (Ss) (Fetter 2001; Tizro et al. 2010). The conventional ways to estimate these parameters are through pumping tests, slug tests, permeameter measurements and grain size analysis (Fetter 2001). However, these tests can be expensive, time-consuming and depend on the existence of boreholes and sampling points conveniently distributed over a study area (Soupios et al. 2007; Perdomo et al. 2018). Likewise, values of hydrogeological parameters can vary over relatively short distances in heterogeneous aquifers, making the results from conventional methods valid only for small sectors of an aquifer (Perdomo et al. 2018). The correlation between resistivity and hydrogeological properties in porous aquifers can be used to get cost-effective estimations of parameters like K and T from surface geoelectrical measurements (Tizro et al. 2010).

The porous aquifer to the north of Oruro (∼300,000 inhabitants) supplies the total demand for domestic consumption, irrigation and industry. However, the hydrogeological knowledge of this reservoir needs to be improved in order to have a future management plan of water resources in the region (D'Elia 2013; GITEC & COBODES Ltda. 2014). The aquifer has been investigated at small scales in the past, aiming to implement wells connected to the supply system (SCIDE et al. 1996). The information from these studies is partially available and corresponds only to a small part of the study area. With the purpose of obtaining preliminary aquifer parameters in an area bigger than where the wells are located, this study combines the available information from drilling reports, pumping tests and geoelectrical measurements to obtain values of K and T, based on an empirical linear petrophysical relationship between resistivity and hydraulic conductivity.

In order to expand K and T estimations, resistivity and aquifer thicknesses were obtained from inverted models corresponding to field measurements applying electrical resistivity tomography (ERT) and transient electromagnetic soundings (TEM), conducted between 2013 and 2017 as part of this study. Complementary data were retrieved from previous studies applying vertical electrical soundings (VES) (e.g., SCIDE et al. 1996; Dames & Moore Norge 2000; Canaviri 2011). The typical resistivity range corresponding to porous aquifers saturated with fresh water was used to estimate hydrogeological parameters, making wider use of the geoelectrical measurements in the study area. This study intends to improve the comprehensive knowledge of one of the most important groundwater reservoirs in the Bolivian Altiplano.

The study area is characterized by a semiarid climate and scarce vegetation. It belongs to the Poopo enclosed catchment where the lake of the same name receives the discharge of surficial and sub-surficial water bodies. However, this lake dried up in 2015 evidencing the decline of available water in the region (Satgé et al. 2017). In terms of groundwater, the most important reservoirs are porous aquifers lying beneath peneplain terrains surrounded by mountains, covering approximately 500 km² at an elevation of ∼3,700 m above sea level (Figure 1). The selected study area is located in the alluvial fan of River Paria. The recharge of the aquifers mainly occurs by lateral infiltration of surface water coming from the rivers, which previously receive and transport the runoff from the mountains during the rainy season. The natural flow direction in the top porous aquifer was mainly to the south, in accordance with the topography. However, the large volumes extracted from a well field in the middle of the alluvial fan have been changing that original flow direction, lowering the water table around the wells and creating a cone of depression of about 5 km radius (Swedish Geological AB 1996; Banks et al. 2002; GITEC & COBODES Ltda. 2014). The study area is considered hydraulically stressed due to the increasing abstraction of groundwater, ∼9.5 million m³/yr (SELA 2018), which in the future might threaten the hydrodynamic equilibrium in the whole aquifer system (D'Elia 2013), since the recharge, ∼12.5 million m³/yr coming from the River Paria, has also been declining (GITEC & COBODES Ltda. 2014).

The geological units in the region consist of three main groups, from older to younger: sedimentary rocks (consolidated), volcanic outcrops (consolidated) and granular deposits (unconsolidated). In the study area, sedimentary rocks including sandstones, siltstones, shales and quartzite are part of the bedrock. The secondary porosity in this unit is assumed the most important for water storage. Overlying the latter, unconsolidated deposits of clay, silt, sand, gravels and pebbles are found in complex arrangements due to the overlapped depositional events (fluvial, lacustrine, glacial and colluvial) (Rigsby et al. 2005). The primary porosity is predominant in these formations, which are considered interconnected as a single porous aquifer (Dames & Moore Norge 2000; D'Elia 2013). A large part of the study area is covered with a thin clay layer which reduces the vertical infiltration (D'Elia 2013).

The flow of water and the flow of electrical current in porous media are governed by similar petrophysical parameters (e.g., porosity, grain size distribution, pore size distribution, mineralogy of grains, geometry of pores, tortuosity), making them related to each other (Soupios et al. 2007; Kirsch & Yaramanci 2009). Many authors have used the analogy of the governing factors for hydraulic and electrical conduction to establish empirical relationships and models (Börner 2009; Perdomo et al. 2018). For example, Archie's law relates bulk resistivity and fluid resistivity in a formation factor that depends on the tortuosity and the porosity (Glover 2016; Perdomo et al. 2018). Another approach is the Kozeny–Carman equation, which relates the hydraulic conductivity among the porosity with its hydraulic radius, the tortuosity and other structural parameters (Börner 2009). In this study, an empirical linear relationship is used to estimate hydraulic conductivity from resistivity, in the spirit of other field investigations that found an inverse correlation between these parameters in different groups of sediments (Mazáč et al. 1985; Kirsch & Yaramanci 2009).

In order to find the relationship between K and ρ valid for the study area, data from previous studies were gathered. A few hydrogeological studies have been conducted around the aquifers near Oruro in the past; however, not all of them are available and the data usually lack precise location, complete record sets and calculations. Information from some boreholes and drilling protocols include stratigraphy, pumping tests and hydrochemical analysis (SCIDE et al. 1996). The mud rotary method was commonly used for drilling. Samples of cuttings are made each metre to get the stratigraphy of the unconsolidated sediments (presented as proportions of gravel, sand and clay). The boreholes reach depths in the range of 60 to 150 m, completed within the porous formations. In order to estimate hydraulic properties, pumping tests were performed at some boreholes and the reported K values are in the range of 5.5 × 10⁻⁵ to 2.9 × 10⁻⁴ m/s (SCIDE et al. 1996; Canaviri 2011). The hydrochemical analyses encompass major ion concentrations and in situ measurements, like pH and electrical conductivity (EC) (Lizarazu et al. 1987; Swedish Geological AB 1996; Gómez et al. 2016). The groundwater salinity in the study area, expressed as EC, is in the range of 0.8 to 3.5 mS/cm; however, only a few sites close to a fault reported the highest values, with the rest of the boreholes reporting values of ≤1.0 mS/cm.

Recent ERT and TEM measurements around the study area achieved good resolution and depth of investigation (DOI), confirming both methods are suitable for the local conditions (e.g., Larsson 2016; Broman & Svensson 2017). Although the methods have different sensitivity degree and resolution (Boiero et al. 2010), both were used to retrieve resistivity values and thickness of the porous aquifer from their inverted models. The ERT method evaluates resistivity by injecting currents and measuring potentials. The measurements using multi-electrode and multi-channel equipment, improve this method covering large distances with good resolution (Dahlin & Zhou 2006). 2-D resistivity sections can be obtained after applying an inverse numerical modelling; an iterative process aiming to reduce the differences between modelled data coming from a synthetic model and the field measurements (Dahlin 2001; Loke et al. 2003). The TEM technique also uses electrical currents, but in this case to generate primary electromagnetic fields that induce eddy currents passing through the ground. The signals of these currents are detected and transformed into resistivity (Christiansen et al. 2009; Reynolds 2011). 1-D resistivity models can be obtained after a numerical modelling process consisting of curve matching to solve inverse scattering problems (Nabighian 1991). The measurements were conducted in different campaigns between 2013 and 2017 using ABEM equipment: the Terrameter LS for ERT and the WalkTEM for TEM.

Finally and to expand the hydraulic conductivity estimations from resistivity, the study includes VES measurements and results reported in some previous studies (e.g., Swedish Geological AB 1996; Dames & Moore Norge 2000; Canaviri 2011; GITEC & COBODES Ltda. 2014). Similarly to ERT, the VES method estimates the resistivity distribution in the ground by injecting currents and measuring potentials, but in this case the models are solved in 1-D. VES data do not follow standard procedures for data acquisition, inversion or results presentation, since the measurements were conducted at different times, using different equipment and inverted with different software. The distribution of the available data is shown in Figure 2. In order to handle consistent 1-D resistivity data from all the available sources, representative points and their correspondent vertical resistivity profiling, were taken from the ERT sections (originally solved in 2-D) (i.e., triangles are used for the further calculations instead of lines marked as ERT).

To estimate hydraulic parameters from geoelectrical tests using a linear relationship, the target aquifer should be clearly visible in terms of resistivity and thickness from the inverted models. For example, in the northwest part of the study area, the ERT profile A is coincident with the resistivity stations 1 and 2 (see Figure 2). As shown in Figure 3, the top part of the section until ∼3,650 masl (∼70 m thickness) is assumed as the porous aquifer with low resistivity signature (∼20 ohm-m) overlying a more resistive layer (∼100 ohm-m) corresponding to the bedrock with low or negligible water content. Two VES measurements are reported at the same sites of stations 1 and 2, their corresponding resistivity models are shown in Figure 3 as columns. These VES models also expose changes in resistivity from low values in the shallow part to higher values at similar levels shown in the ERT model. Despite this change in resistivity supports the inferred change in geological facies, the values shown in the VES models are not necessarily the same as on the ERT model due to the differences between these methods (different equipment, sensitivity, processing and inversion procedures). In this case, the resistivity coming from the ERT is assumed more reliable since the raw data were processed in this study, whereas the VES models were obtained from a thesis showing model results directly without raw data (Canaviri 2011), hence, the data quality of that study could not be assessed. The resistivity and target aquifer thickness for stations 1 and 2 come from the ERT line A.

The salinity of water influence greatly the resistivity measurements regardless the method applied to get that information. The top porous aquifer yields fresh water (EC <1.0 mS/cm), but the fractured aquifer might be saturated with saline water around some faults (EC >3.5 mS/cm). The knowledge about the fractured aquifer is limited; however, its high salinity and temperature can be attributed to the existence of deep thermal sources, most likely the volcanic intrusions mentioned in the regional geological background. Not all the faults in the region seem to reach thermal sources, but those that presumably do, have a distinctive low resistivity signature. The hydrothermal flows seem to go upwards in some parts along the faults; consequently, some sectors of the porous aquifer present hydrothermal intrusions with artesian conditions. For example, one of these faults with hydrothermal flows crosses the study area coming from the south. The well WP10 (see Figure 2) is apparently located along the projection of the mapped fault. The drilling depth in this well is ∼100 m and the water is more saline (EC ∼3.5 mS/cm) than in others in the vicinity. Although the salinity in WP10 is greater than the limit for consumption (1.5 mS/cm) (IBNORCA 2010), mixing this source with others from the porous aquifer keeps this parameter below the limit before being distributed. The closest ERT resistivity section (B) shows a feature with low resistivity that is assumed as corresponding to the hydrothermal intrusion (Figure 4).

The DOI of the ERT model B is ∼220 m below the surface. The section seems completed within the porous aquifer, in a site that is also indicated as the deepest part of the porous aquifer according to previous studies (e.g., Banks et al. 2002). However, the hydrothermal intrusion precludes the use of the linear relationship due to the saline intrusion. In addition, the very shallow layer of the section has relatively high resistivity values (∼80 ohm-m); it is interpreted as the part above the water table, i.e., unsaturated sediments. The same feature also appears on the ERT line C, but its resistivity value is slightly lower (Figure 5). Unlike ERT sections B and C, the unsaturated fringe does not appear in section A presumably because it is located outside of the influence of the cone of depression shown in Figure 1; therefore, the water table around section A might be shallow.

The target aquifer is also visible on the TEM inverted models. For example, station 28 (Figure 6) has the three main parts of the hydrogeological settings in the study area: the unsaturated fringe on top, the porous aquifer and the bedrock, each one with a distinctive resistivity.

The empirical relationship between K and ρ can be established using the hydraulic conductivity obtained from pumping tests conducted in wells located in the alluvial fan of River Paria, and found in technical reports (e.g., SCIDE et al. 1996; Dames & Moore Norge 2000). The corresponding resistivity was adopted from the closest geoelectrical station to each well (see Figure 2). Figure 7 shows K vs. ρ values fit in a linear correlation except for WP10. As mentioned before, the water in well WP10 is influenced by high salinity associated with hydrothermal intrusions because it is located along the projection of a fault that presumably reach thermal sources. The salinity of water in this place increases the electrical conductivity of the ground measured with VES (Canaviri 2011), making this point out of the linear tendency. The rest of the points were used to establish a linear relationship indicated in Figure 7, which serves to predict K values in the rest of the study area based on ρ values from geoelectrical models.

Kirsch & Yaramanci (2009) present a similar linear relationship for a fine grained sand aquifer, which type of sediment is likely abundant in the present study area. Hydraulic conductivity in the rest of the geoelectrical stations were estimated using this linear equation. Furthermore, transmissivity was also estimated using the relationship T = K · h (Fetter 2001), where h is the thickness of the saturated porous aquifer; in some sites it is underlying an unsaturated fringe (as shown in Figures 4–6). The saturated thickness as well as resistivity were obtained from inverted resistivity models; in the case of VES and TEM these models are solved in 1-D, and ρ and h values can be directly obtained from them (as in Figure 6). For ERT models solved in 2-D, ρ and h values were obtained vertically (like 1-D) below the location of each station from the blocked inverted models, regardless of lateral changes in resistivity in the rest of the section. The last column T_reference is obtained multiplying K_reference by h.

The distribution of estimated T over the study area allows the most promising areas in terms of groundwater production to be identified; these values were interpolated with the Kriging method to obtain smooth transitions, as shown in Figure 8. The highest transmissivity is located in the middle of the study area, in agreement with previous studies indicating the same sector holds the thickest portion of the porous aquifer (e.g., Banks et al. 2002). On the other hand, the lowest values are located in the southwest part, where the top of the bedrock was detected at shallower levels and where the hydrothermal intrusion constrains the usable thickness around the fault shown in Figure 2. Interpolated values in the northern top and southeast parts of the study area may not be entirely reliable due to the less densified resistivity stations for interpolations.

The estimation of aquifer characteristics presented in this study is based on the relationships between hydraulic and electrical properties in unconsolidated sediments. The linear correlation found between K_reference (obtained from pumping tests) and ρ (obtained from inverted geoelectrical models) is valid only for the study area, since this type of inverse correlation has been empirically established varying according to the rock type and the grain size (Mazáč et al. 1985; Kirsch & Yaramanci 2009). The factors governing the geoelectrical response are also critical when it comes to evaluating the validity of the method. In this aspect, the biggest uncertainty comes from the mixed use of results from ERT, TEM and VES, since these methods evaluate resistivity in different ways and have different sensitivities. However, in the spirit of making wider use of the available resources in the study area the proposed empirical petrophysical relationship can provide results worth analysing.

Regarding the thickness of the porous aquifer, in some parts of the study area the contact between unconsolidated sediments and bedrock can be identified on the resistivity models as transitions from relatively low values (∼20 ohm-m) to greater values (e.g., ERT section A, Figure 3). Likewise, saline water intrusions can be seen as transitions to lower values (e.g., ERT section B, Figure 4), reducing the thickness of the aquifer saturated with fresh water. In cases where the bottom limit of the target aquifer cannot be identified due to limitations in the DOI of the geophysical methods, the results might have a degree of uncertainty. In the central part of the study area, stations 7 and 14 were assigned the maximum thickness of their respective DOI (see Table 1), although the thickness in these points may be greater in reality. Nevertheless, the same site was identified as having the deepest unconsolidated sediment sector in previous studies (e.g., Dames & Moore Norge 2000; Banks et al. 2002).

The density of geoelectrical stations plays a significant role in resolving hydrogeological properties in new areas. The greater the distance between stations, the less reliable is the interpolation of intervening regions. The access to some parts of the study area was restricted during the fieldwork due to conflicts between villagers and water well owners. Consequently, it was not possible to extend the measurements to cover blank spaces in between geoelectrical stations shown in Figure 8. The estimated transmissivity values in the study area vary between 3.4 × 10⁻³ and 4.0 × 10⁻² m²/s evidencing, in turn, the influence of the great variations in the aquifer thickness.

The central part of the alluvial fan of the River Paria is favourable for geoelectrical surveys given its clear distinction between geological facies, greatly influenced by the type of water saturating the aquifers. The top porous aquifer, used for fresh water supply, was evaluated in its hydraulic properties by using an empirical relationship between resistivity and hydraulic conductivity to identify suitable sites for groundwater production, which, in this case, is the central part of the investigated area.

Resistivity models from different geoelectrical measurements (electrical resistivity tomography, transient electromagnetic soundings and vertical electrical soundings) were evaluated in terms of target aquifer thickness and resistivity from the inverted models. The distribution of the calculated transmissivity in the study area is also in good agreement with geological models from previous investigations; the highest values are located in the centre of the study area, at the same site where the porous aquifer is thicker. In sites where the top of the bedrock is shallow and where hydrothermal intrusions were detected, the estimated transmissivity has lower values.

Getting an early idea of hydrogeological properties like transmissivity and hydraulic conductivity, before conducting more demanding and costly tests, can be very valuable, identifying the most promising sectors for groundwater production. A more ambitious study may include extended surveys to cover larger areas (e.g., with airborne geophysics), necessary for groundwater management plans.

The study was funded by the Swedish International Development Cooperation Agency (SIDA) and the Society of Exploration Geophysicists (SEG) through the Geoscientists Without Borders (GWB) programme. This work was also supported by Lund University (LU) in Sweden, Aarhus University (AU) in Denmark, Universidad Mayor de San Andrés (UMSA) and Universidad Técnica de Oruro (UTO) in Bolivia. The TEM equipment was provided by Engineering Geology Division LU with support from Guideline Geo AB. The ERT equipment was provided by UTO. We are grateful to the people who helped to conduct the field measurements: Ramiro Pilco, Måns Larsson, Emil Svensson, Kamran Ahmad, Mauricio García, Rafael Mendoza, personnel of Instituto de hidráulica e hidrología (IHH-UMSA), Servicio Local de Acueductos y Alcantarillado de Oruro (SELA), Unidad de Saneamiento Basico y Vivienda (UNASVBI), Ministerio de Medio Ambiente y Agua (MMAyA) and Corimex Ltda., students from UTO and many other volunteers. Finally, we acknowledge the assistance of the facilitators at the workshop on scientific writing, arranged by Aalto University, Helsinki, Finland and IWA/YWP Finland in August 2018. This has considerably improved the content of the paper.