Finite Element Analysis and Experimental Characterization of Soil Electrical Resistivity at El Roque de los Muchachos Observatory ()
1. Introduction
El Roque de los Muchachos Observatory (ORM) is located at the rim of Taburiente National Park, at heights ranging from 2340 to 2396 meters above sea level, in the municipality of Garafía (La Palma, Spain). It hosts one of the largest fleets of telescopes in the world, including infrared, visible, microwave and gamma-ray instruments. ORM sky quality continuously attracts projects aimed at building large scale instruments. The Cherenkov Telescope Array (CTA) project is currently the most ambitious one, with plans for building four Cherenkov Large Size Telescopes with a 23 m dish diameter and a subarray of Medium Size Telescopes with 12 m dish diameter [1] [2]. The construction requirements approved for these instruments include tough specifications on the quality of the earth resistance, for which the terrain resistivity is a crucial factor. These specifications are mainly motivated by the need to comply with safety regulations as well as to mitigate the risks of damaging sensitive instrumentation due to static discharges and lightning strikes. However, no information about the electrical resistivity of the soil at ORM was found.
The four-point Wenner method [3] is currently the most popular procedure for electrical resistivity measurements, due to its simplicity and the low cost of the required equipment. Wenner’s original work provides a simple equation to derive the soil resistivity ρ from the measured current I and voltage V, and the interelectrode distance a in a four-point configuration,
. (1)
This expression holds provided that several conditions are fulfilled, namely:
1) The resistivity is homogeneous;
2) The electrodes are point-shaped and placed on the surface, being interelectrode distance much larger than interelectrode depth;
3) The conductive volume is semi-infinite;
These measurements can only give what is usually called apparent resistivity, since in most cases the result is obtained from the contribution of different local resistivities within the volume through which the current flows. Despite these limitations it has been shown that Wenner method still provides accurate results even if some of the requirements for the model validity are not fully satisfied [4].
On the other hand, Schlumberger method provides a convenient alternative to speed up the measurements when detailed profiling is needed, since the resistivity can be obtained in several points by moving only the two electrodes of the voltage channel [5]. It is also convenient when using instrumentation lacking accurate voltage measurements and large interelectrode distances are needed. This is due to the fact that the Schlumberger method allows you to place the voltage electrodes closer to the current ones than in Wenner method [6]. This way the voltage to be measured is higher. Special procedures have also been proposed to determine soil resistivities in thin layers [7].
Other Geological surveys regularly use techniques like Electrical Resistivity Tomography (ERT) for hydrogeological, mining and geotechnical investigations [8] [9]. These techniques use electrodes arrays with the same basic operating principles as the Wenner and Schlumberger methods. The patterns are built by means of numerical, inverse problem algorithms. 3D tomography techniques are especially powerful to characterize low resistivity (typically below 10 kΩ·m) soils when the terrain allows the proper insertion of a large number of electrodes in specific locations.
Mathematical expressions for the apparent resistivity in earth structures with an arbitrary number of layers (single, double or multilayered) have been derived for Wenner method [10] [11] [12]. Finite Element Modeling (FEM) software has been used in previous works to study the influence on electrical resistivity measurements in materials. However, most of these studies are focused on the role of concrete presence on soil resistivity [4] [13] [14]. FEM has been demonstrated to be an alternative method to realistically model the boundary conditions and the electrode geometry versus analytic methods, which rely on unnecessary approximations.
The original Wenner method was found to be the most convenient procedure to characterize experimentally the electrical resistivity of the ORM soil due to the optimum balance it offers between the degree of detail obtained and the required cost and time efforts. Figure 1 shows two images of the measured areas.
The ORM is in a protected area where the soil is largely covered by adenocarpus viscosus plants, an endemic species which can reach heights of up to 1.5 m with deep and hard roots. In addition, the soil of this site is mostly filled of hard rocks and materials with resistivities exceeding 10 kΩ·m in some locations. Other measurement methods, such as 3D tomography, were ruled out due to the difficulty in carrying them out.
In this work, we present measurements done in six different locations of an area at ORM covering around 250 m × 275 m and centered at coordinates 28˚45'42.9"N, 17˚53'28.5"W. These locations correspond to the sites where the telescopes MAGIC I, MAGIC II and LST1 are built, as well as the future locations of LST2, LST3 and LST4. The resistivities at the sites have been simulated with COMSOL Multiphysics software using two different models: a simple single layer model and a three-layer model. Practical limitations of the measurements and the consistency of the mathematical model with the true terrain structure are discussed with the support of previous geotechnical studies and available literature of electrical properties of the materials of which the soil is composed.
![]()
![]()
Figure 1. Views of the LST area landscape showing the dominant vegetation.
2. Experimental Procedure
2.1. Resistivity Measurements
The resistivity meter used for the measurements is the Chauvin Arnaux’s C.A 6470N Terca 3 [15], with its corresponding 150 m earth and resistivity kit’s pikes and wires. This equipment, among other applications, calculates the apparent resistivity using Wenner method.
Wenner method provides a value of the average apparent resistivity at a point located around the middle of the pike line and at a depth of the order of the pike distance [16]. The standard Wenner protocol with occasional support of special procedures to improve the contacts, like pouring small amounts of saltwater around the pikes, worked properly at all sites except at LST3 area. The measurements at this area were particularly difficult to obtain, and several trials were needed to find places for the pikes giving reliable results. However, the measurements became significantly easier when the pike separation was increased. In order to understand this behaviour, it is necessary to revise the instrument capabilities. Table 1 shows the relevant specifications of the equipment [15]. As it can be observed, the maximum resistivity the instrument can measure is around 1 MΩ·m, which is well above the maximum reliable values measured onsite at all distances. The test signal was 32 V and the AC frequency 128 Hz.
Table 2 provides the threshold values which trigger an error message in the instrument. The one which triggered the observed errors in the display during
![]()
Table 1. Characteristics and limitations of the equipment [15].
![]()
Table 2. Error indications according to manufacturer. 3P and 4P refer to instrument functions that were not used for the results shown in this work. IH-E is the measured current applied to current path. Ui correspond to effective measured voltages between the electrodes.
the short-range measurements at LST3 is shown in the first row, which indicates that the current measured in the current channel IH-E was too low. The lack of other errors ensures us that stable contacts were achieved, and the measurements were free from parasitic interferences and stray signals.
We therefore associate these errors to the own irregularities of the terrain. We know that it is composed of different layers, but we believe that it is possible that within these layers there are also irregularities that cause the terrain to behave discontinuously, thus worsening the circulation of the current IH-E. This would prevent the measurement from being taken correctly, as the current path inside the ground does not drive as it should. This behavior is illustrated in Figure 2, where we see the current represented as a dotted line, and a possible obstacle in the circulation of it and, therefore, in the closing of the circuit. Only the separation of the pikes could enable the generation of current paths below the obstacle and therefore the measurement of detectable currents. All unreliable measurements were disregarded for the modeling presented here.
For this study, two measurement campaigns have been carried out, in September 2018 and May 2019 respectively, from which the data extracted are displayed in Figure 3. MAGIC I, II and LST1 measurements were made around a terrain that was altered for the installation of the telescope foundations. The construction works required the removal of land for the foundation of each telescope
![]()
Figure 2. Schematic representation of a four-probe measurement with a current cut produced by a non-conductive obstacle. The obstacle is represented by the grey rectangle below blue pike, and current is represented by a dashed line.
![]()
Figure 3. Soil resistivity measurements made in ORM. The numbers correspond to the interelectrode distance in m.
up to a depth around 2 m. The humidity values in the area were monitored during the measurements with the weather station facility of the MAGIC telescopes. Recorded relative humidity remained in the range [30% - 41%] for all measurements. At least three measurements with different pike distances were used at every site in order to check if the resistivity varies with depth.
Average resistivities were 1.44 kΩ·m, 5.14 kΩ·m, 2.45 kΩ·m, 0.57 kΩ·m and 0.53 kΩ·m for LST1, LST2, LST4, MAGIC I and MAGIC II, respectively. LST3 site exhibited resistivities ranging from 1.74 kΩ·m to 11.5 kΩ·m. The fact that resistivity decreases as interelectrode distances increase in this site indicates that there is a high resistance layer in the upper part of the soil whereas a layer of a lower resistivity was found at depths around 20 m.
2.2. Soil Composition
Additional information on the terrain can be obtained from Geotechnical Studies (GTS), where the physical evaluation of the soil is made, to identify its composition. So, these excavations give data about the subsurface material at different depths. The soil where the telescopes will be located in had to be analysed according to the “Guide for planning and carrying out Geotechnical Studies for building in the Comunidad Autonoma de Canarias” (GETCAN-011 [17] ).
It is known that La Palma, part of Canary Islands, has a volcanic origin as almost the 100% of the geology of this archipelago. In particular, for the LST1-LST4 area we have altered basalt masses. The geotechnical unit is n˚ III, based on public geological and cartographic databases [18] [19], which explain that these types of soils correspond to basic castings of small thickness and moderate high alteration, which manifest as vertical alternation of basaltic compact levels.
Like resistivity measurements anticipate, the soil composition obtained in GTS where LST1 is placed, is almost a homogeneous layer of slag altered with clay (≈15%). However, the LST3 soil has layers of several materials. A reconstruction of the terrain based on the excavations made in different points around LST3 site is shown in Figure 4.
The excavation in point T2 is the closest to the resistivity measurements site in LST3 locations. According to this, composition of the soil in LST1 and LST3 is shown in Table 3. The deepest layer of slag altered with silts is not considered since it was found only in T2 survey. We also neglected the deepest layer of fractured basalt for the same reason.
We can conclude that, in general, the rocks of LSTs area are a mixture of igneous and volcanic types [17]. According to the composition of the rocks, they are MAFIC with other materials (clay, gravel and silts) [20].
Knowing the composition of the soil, one would expect that the representative resistivity values can be approximated based on the electrical properties of each rock type. Therefore, the resistivity of the rock types found in the GST of LST1 and LST3 areas have been checked from several references and collected in Table 4.
![]()
Figure 4. Draft of LST3 study based on the excavations made in different points (T1, T2 and T3) around the telescope.
![]()
Table 3. Composition of the soil based on the GTS closest to the measurement sites (T2 location).
![]()
Table 4. Representative resistivity values of the rock types found in GTS.
It can be observed that there is a considerable difference among the resistivities reported by references [21] - [26] for the materials found in the ORM GTS. References [21] and [25] provide a narrower margin of resistivities for each material than the others because they use a more detailed classification. Data shown in Table 4 from these references correspond to those types which are closer in nature to the ORM soil. On the other hand, resistivities reported in reference [25] offer the advantage that they distinguish between basalt and slag or scoria. There is a wider clasification of the sedimentary rocks too. Reference [23] considers all possible igneous rocks as a single group.
It can be concluded that no reliable data on resistivities could be found in available literature. Furthermore, the composition of each soil layer is not necessarily homogeneous, and other factors like the humidity or cracks have a deep impact on the resistivity. Therefore a dedicated resistivity survey onsite is essential.
3. Finite Element Model
FEM has been used to calculate the apparent resistivity equipment using Wenner method should measure in two representative locations, which were selected based on the ORM measurements. According to the results shown in Figure 3, five sites (MAGIC I, MAGIC II, LST1, LST2 and LST4) exhibit a homogenous resistivity. LST3 exhibited a different behaviour, with a strong dependence of measured resistivity with pike separation. Therefore, two different simulation models have been implemented. We define a “Soil 1” type as a homogeneous soil. “Soil 2” type corresponds to a multi-layered soil, as the one found in LST3. The aim is to find by simulation of the electric field and current density the true terrain resistivities which would reproduce the measurements displayed by the equipment when the interelectrode separation distance changes in the range of interest, from 2 m to 30 m.
The model uses COMSOL Multiphysics v5.3a software. For this case, Electric Currents Module and a parametric sweep for a parameter in the frequency domain have been used for the study. The computational domain is a 3D parallelepiped containing the complete geometry modeled.
The electrodes are considered as “copper material” defined in COMSOL. Material parameters needed for each soil layer are displayed in Table 5. These data are based on the GTS made in LST1 (Soil 1) and LST3 (Soil 2). Several COMSOL studies have been made to estimate the resistivity values for each layer which lead to the best agreement with the measurements.
To avoid contour artifacts due to the finiteness of the model volume, the dimensions of the parallelepiped have been chosen large enough: its length is seven times the maximum electrode separation, 30 m, and its width is 3.5 times this distance. The parallelepiped has 65 m depth. An Infinite Element Domain condition has been imposed to an external layer of 2 m thickness in
![]()
Table 5. Soil resistivity parameters for the multi-layered soil designed in the COMSOL models.
order to improve the simulation of an infinite volume.
The electrodes were simulated as square patches of 10 cm side. Although they could be simulated as point sources as well, the details of their geometry are irrelevant due to their negligible size as compared to interelectrode spacing [4]. For boundary conditions, current path electrodes are considered as voltage sources of ±16 V, respectively, which correspond to the 32 V source applied during the measurements, as indicated in Table 1. Ground is located at the limits of the parallelepiped. The simulation provides the voltage difference V between the inner electrodes and the current flowing thru the outer ones I. Then, we calculate the apparent resistivity applying Equation (1).
In our models, the mesh consisted of tetrahedral elements with a nonhomogeneous density. The extra fine default condition set by COMSOL was chosen for all the volume except a local region around the electrodes. A cube of 2 m side around each electrode was set to extremely fine mesh density condition. Figure 5 shows the details of the implemented model geometry.
4. Results and Discussion
The results of the simulation are shown in Figure 6 for a soil type 1, represented by the LST1 site, and Figure 7 for a soil type 2, corresponding to LST3. Each figure displays an image with an aerial view showing the place of the GTS survey as well as the line along which resistivity measurements with stable readings could be made. Also shown are the volume geometry simulated, a plot with both equipotential surfaces and current lines and the comparison between simulated and measured values of the apparent resistivity.
In the resistivity plots we include three simulations for each case, in order to provide an insight of the influence of the first layer characteristics on the final result. In soil 1 case (LST1 site) it is found that the best fit is obtained for a resistivity of 1.44 kΩ∙m. Simulations agree with measurements within an error of
![]()
Figure 5. Implemented model geometry. Soil (multilayered case) (1). Infinite Element Domain Layer (2). Extremely fine mesh zone surrounding electrode (3). Electrodes patches (4). Ground (5).
4.1%. In soil type 2 the thickness of the first layer h makes a strong influence on the final result. For the adequate thickness (h = 2.85 m), discrepancies with measurements are close to 10%.
It can be observed how the current density curves between electrodes bend if the terrain is not homogeneous, as also happens with the equipotential surfaces, which lose their shape of perfect hemispheres. Therefore, the ideal behavior that is often assumed in theory does not work for these cases.
The measurements made at LST3 for the closest distances provided values around 12 kΩ·m, which were well below the maximum resistivity limit of the instrument. Nevertheless, there were several cases in which the measurement could not be done. This led us to the conclusion that the terrain at LST3 exhibits shallow obstacles providing current cuts that make short range readings difficult or impossible to obtain. The nature of these obstacles could not be clarified with the available GTS surveys, but despite the difficulties a layer of low resistivity, around 1.5 kΩ·m, was found when measurements with large interelectrode distances were made. So, if electrode separation is increased and the measurement corresponds to depths where the conductivity is greater, current lines flow through this area and the reading becomes feasible.
According to the simulations, we can also confirm that the main responsible for the low resistivity found at high depths could be the layer of volcanic agglomerates of slag and gravel located at depths around 8 m and below. The value of soil resistivity is determined by the properties of the terrain components: rock type or nature, the amount of salts dissolved, the water content, the compactness… Moreover, available literature on electrical resistivities of similar materials exhibited large discrepancies among the authors and could not provide us reliable data. Nevertheless, the orders of magnitude of resistivities found in this work are compatible with the wide margins reported in references [22] [23] [24] and [26], although they are not in agreement with the data provided in [21] and [25].
5. Conclusions
This work has presented a study of the electrical behavior of the terrain at El Roque de los Muchachos Observatory, which covers six different places located within an area centered at coordinates 28˚45'42.9"N, 17˚53'28.5"W and with an extension around 250 m × 275 m. This study has revealed a significant variability of the terrain not only in terms of the average apparent resistivity but also in the own layer structures.
While five of the studied places showed a behavior compatible with a single layer structure up to depths of around 10 m, the sixth one, the LST3 site, exhibited a dramatically different behavior. The soil at LST3 contained materials in the shallow layers with significantly higher resistivities, and it was considerably difficult to find places for the pikes that provided stable readings.
FEM simulations could provide an insight on the correlation between materials found in the GTS surveys and the measured resistivities. The four-layer model used for LST3 was found to accurately predict the measured resistivities while keeping a reasonable compatibility with the available surveys. The simulations confirm resistivity values ranging between 15 and 20 kΩ·m for the fractured basalt layers found at El Roque, with lower values for deeper layers (<2 kΩ·m).
We believe that this information might help to efficiently optimize the future ground structures for the telescopes under construction at ORM or other similar sites with soils of volcanic nature.
Acknowledgements
This work has been funded by the Spanish Ministry of Science and Innovation, thru the research project FPA2017-82729-C6-4-R.