Hydrogeophysical Assessment of the Critical Zone below a Golf Course Irrigated with Reclaimed Water Close to Volcanic Caldera

: The geometry and the hydraulic properties of the unsaturated zone is often difﬁcult to evaluate from traditional soil sampling techniques. Soil samples typically provide only data of the upper layers and boreholes are expensive and only provide spotted information. Non-destructive geophysical methods and among them, electrical resistivity tomography can be applied in complex geological environments such as volcanic areas, where lavas and unconsolidated pyroclastic deposits dominate. They have a wide variability of hydraulic properties due to textural characteristics and modiﬁcation processes suh as compaction, fracturation and weathering. To characterize the subsurface geology below the golf course of Bandama (Gran Canaria) a detailed electrical resistivity tomography survey has been conducted. This technique allowed us to deﬁne the geometry of the geological formations because of their high electrical resistivity contrasts. Subsequently, undisturbed soil and pyroclastic deposits samples were taken in representative outcrops for quantifying the hydraulic conductivity in the laboratory where the parametric electrical resistivity was measured in the ﬁeld. A statistical correlation between the two variables has been obtained and a 3D model transit time of water inﬁltration through the vadose zone has been built to assess the vulnerability of the aquifers located below the golf course irrigated with reclaimed water.


Introduction
Golf courses irrigation using reclaimed water provides a significant and viable opportunity to ensure the supply, sustainability and resilience of local water resources [1,2]. There is an enormous potential for treated wastewater use for agricultural irrigation purposes [3,4] but some barriers exist to widespread adoption due to some potential contaminants that have side effects on the earth's critical zone affecting aquifers, the quality of soil, and/or public health [5,6].
Generally, precise information about the spatial variation at field-scale soil hydraulic properties is essential to carry out a careful exploration of the critical zone [7]. The subsurface geology guides the water movement specially after large rainfall events. As these events occur frequently in arid and semiarid zones, subsurface knowledge is a critical factor to determine water management guidelines. Traditional hydrological methods are

Study Area
The hydrogeophysical study was carried out in a golf course located midlands of Gran Canaria island's north-eastern section, at an altitude of between 400 m and 500 m ( Figure 1). The Bandama Golf Course has 18 holes whose fairways and greens cover approximately 14.5 ha and spraying irrigation frequencies vary between winter and summer when doses reach a maximum rate of 7 mm/day [19]. From the climatological point of view, the Bandama Golf Course is in an area with an annual rainfall module slightly above 300 mm, while the average temperature is 19 • C (22 • C in summer and 16 • C in winter).
The rocks that outcrop in the area are Holocene basaltic lava and pyroclasts. These materials were emitted in the most recent eruption of Gran Canaria (1970 ± 70 Before Present), where a strombolian cone (Pico Bandama) and a phreatomagmatic caldera (Caldera de Bandama) arised. Pyroclastic deposits consisted of tephra air fall deposits and pyroclastic flows (surges) covering a surface of 50 km 2 [22,23]. The Caldera of Bandama, is 900 m in diameter and 250 m deep, and the golf course is located within its western sector (Figure 1a). Therefore, as Figure 1b shows, the eastern slope of the Caldera allows direct access to the geology of the unsaturated zone composed by: (1) Miocene phonolithic basement that includes interbedded alluvial conglomerates of the Las Palmas Detritic Formation, (2) Pliocene fractured basanitic lava flows and landslide breccia from the Roque Nublo Group and (3) Holocene pyroclastic deposits emitted in the phreatomagmatic eruption of the caldera itself.
Water 2021, 13, x FOR PEER REVIEW 3 of 15 The rocks that outcrop in the area are Holocene basaltic lava and pyroclasts. These materials were emitted in the most recent eruption of Gran Canaria (1970 ± 70 Before Present), where a strombolian cone (Pico Bandama) and a phreatomagmatic caldera (Caldera de Bandama) arised. Pyroclastic deposits consisted of tephra air fall deposits and pyroclastic flows (surges) covering a surface of 50 km 2 [22,23]. The Caldera of Bandama, is 900 m in diameter and 250 m deep, and the golf course is located within its western sector ( Figure 1a). Therefore, as Figure 1b shows, the eastern slope of the Caldera allows direct access to the geology of the unsaturated zone composed by: (1) Miocene phonolithic basement that includes interbedded alluvial conglomerates of the Las Palmas Detritic Formation, (2) Pliocene fractured basanitic lava flows and landslide breccia from the Roque Nublo Group and (3) Holocene pyroclastic deposits emitted in the phreatomagmatic eruption of the caldera itself. Two main different soil types have been characterized depending on their origin. In situ soil is a Torriarents (adjacent natural soils are vitriotorrands) and transported soil corresponds to an Ustalfs dominated zone [24]. The in-situ soil consists of slightly altered basaltic pyroclasts with a thickness of 0.25-0.5 m, on which a sandy-loam alteration cover has developed.
The soil transported from agricultural lands of higher elevations of the same slope of the island was used for the construction of four fairways of the golf course and is identified as silty-clay nature. Recent studies [21] identified a different behavior of both soils and have shown that variability of soil parameters are influenced by irrigation management, soil type, water quality and quantity, and seasonality of sampling. Two main different soil types have been characterized depending on their origin. In situ soil is a Torriarents (adjacent natural soils are vitriotorrands) and transported soil corresponds to an Ustalfs dominated zone [24]. The in-situ soil consists of slightly altered basaltic pyroclasts with a thickness of 0.25-0.5 m, on which a sandy-loam alteration cover has developed.
The soil transported from agricultural lands of higher elevations of the same slope of the island was used for the construction of four fairways of the golf course and is identified as silty-clay nature. Recent studies [21] identified a different behavior of both soils and have shown that variability of soil parameters are influenced by irrigation management, soil type, water quality and quantity, and seasonality of sampling.
Since 2002, the installation of a tertiary desalination treatment system has significantly reduced the salinity of the reclaimed water (1000 µS/cm) and since December 2009, the quality has further improved to 300 µS/cm. This change in irrigation water quality had a direct effect on the parameters measured in the soil and in the water collected in the lysimeters installed in the field, pointing to the destabilization of soil aggregates [21].
The island hydrogeological conceptual model can be sketched as a unique groundwater body recharged by rainfall infiltration and discharged into the sea or some discharging points into springs and ravines. In the area, the aquifer system mainly exploits phonolitic materials using 2.5-3 m in diameter wells reaching depths in the 15-300 m range [10]. The water table is located 250 m below the Bandama Golf Course and groundwater flow from the golf course to the Las Goteras Ravine has been previously defined (Figure 1a). The setting-up of a monitoring network of water points along the ravine has made it possible to characterize the groundwater quality and, also, the presence of emerging contaminants and priority substances in the aquifer [6].

Electrical Resistivity Tomography
An electrical resistivity tomography (ERT) survey was conducted to assess the subsoil properties of the golf course. The method is based on measuring the potentials between one electrode pair while transmitting DC between another electrode pair (quadrupole). The depth range increases with increasing space between the current electrodes, whereas a shorter separation increases resolution [25]. The ERT uses fixed multiple electrodes in the soil surface that change function automatically according to the acquisition array previously selected. All possible combinations of quadrupoles are considered, resulting in a dataset of apparent resistivities at the so-called pseudo-depth at different locations. The large volume of data gathered by multielectrode systems requires automated data handling and processing [26].
ERT data was acquired with a Syscal Pro resistivity meter (IRIS instruments, Orléans, France). The system features an internal 250 W power source and an internal switching board for 48 electrodes. The quadrupole array chosen was Wenner-Schulmberger because it is sensitive to both vertical and horizontal structures and has an adequate signal strength [27]. The array has high performance and stability in high electrical resistivity environments such as volcanic rocks and it is effective for the characterization of horizontal or slightly inclined layers that have lateral facies variations and/or verticalized structures, as is the case of the studied setting [28,29].
RES2DInv was the software used for the inversion of the ERT data and to estimate the true resistivity of the subsoil [30]. The subsurface is divided into fixed dimensions cells and the procedure is based on the smoothness-constrained least-squares method. The resistivity values are adjusted iteratively until a suitable agreement between the raw data and the model responses is reached, based on a nonlinear optimization technique by least-squares fitting [31]. During the inversion procedure, the root-mean-square value of the difference between experimental data and the updated model response is used as a convergence criterion.
In the present paper, the robust method was selected. The method assumes that the subsurface consists of limited homogeneous regions with a sharp boundary among them. The robust scheme is the reasonable choice where the subsurface comprises units with sharp limits to accurately define both layer boundary locations and layer resistivities. Indeed, it produces models by minimizing the absolute value of data misfit, becoming more efficient in removing noise compared to other inversion methods [32].
The design of geophysical surveys has the objective to cover the study area with a representative grid of the variability of electrical resistivity values. The profiles were disposed as regularly as possible in the site and their location was conditioned by the morphology of the fairways and for not disturbing the development of the golf game during the acquisition procedure ( Figure 2).  As a result, we use 48 electrodes arrays to obtain 2D ERT cross-sections with 94 m length, reaching an investigation depth close to 20 m and a resolution of two meters apart between geoelectrical values. The data collection includes 941 quadrupoles for each profile and rs check resistance between adjacent electrodes always below 10kOhm. To validate each measurement, we have repeated it, or stacked, from three to five times requesting a standard deviation for the group of stacked measurements of 3% maximum.
Geoelectrical data was positioning with a differential GRS1 GPS (Topcon, Itabashi, Japan), and relative relief profiles of the cross-sections were converted into georeferenced elevation profiles using an earth digital elevation model provided by the Spanish Geographical Survey (IGN). The elevation model has a 2 × 2 m resolution and the absolute vertical accuracy corresponds to an average mean quadratic error of 0.15 m in flat and low vegetation areas.
The subsequent subsurface characterization must consider the overlapping of resistivity values for different rocks and soils because the resistivity depends on several factors, such as mineralogy, soil water content, grain size distribution and porosity. For instance, clayey soil normally has lower resistivity than sandy soil and an air-filled porosity soil type will have higher resistivity values conversely to a water-filled porosity soil type As a result, we use 48 electrodes arrays to obtain 2D ERT cross-sections with 94 m length, reaching an investigation depth close to 20 m and a resolution of two meters apart between geoelectrical values. The data collection includes 941 quadrupoles for each profile and rs check resistance between adjacent electrodes always below 10kOhm. To validate each measurement, we have repeated it, or stacked, from three to five times requesting a standard deviation for the group of stacked measurements of 3% maximum.
Geoelectrical data was positioning with a differential GRS1 GPS (Topcon, Itabashi, Japan), and relative relief profiles of the cross-sections were converted into georeferenced elevation profiles using an earth digital elevation model provided by the Spanish Geographical Survey (IGN). The elevation model has a 2 × 2 m resolution and the absolute vertical accuracy corresponds to an average mean quadratic error of 0.15 m in flat and low vegetation areas.
The subsequent subsurface characterization must consider the overlapping of resistivity values for different rocks and soils because the resistivity depends on several factors, such as mineralogy, soil water content, grain size distribution and porosity. For instance, clayey soil normally has lower resistivity than sandy soil and an air-filled porosity soil type will have higher resistivity values conversely to a water-filled porosity soil type and it has been necessary to incorporate soil and geological setting to improve the interpretation of the ERT results [33].

Hydraulic Conductivity
Hydraulic conductivity is the key factor of water flow through the substrate and it is affected by in-situ structure and pore volume [34]. Particularly, saturated hydraulic conductivity (Ks) is used to describe the movement of water through saturated soils and is a critical component in a resource management decision such as water conservation and irrigation systems [35]. Saturated hydraulic conductivity has been measured from undisturbed representative soils samples and volcanic deposits taken directly from selected outcrops. More than twenty soil samples from two described top profiles (both from lane and rough) were analyzed in each of the two sampling periods. The collection was carried out by driving a standardized 250 cm 3 cylindrical sampler into the soil. Once in the laboratory, the previously prepared soil cylinders were watered from the bottom until saturation is reached and then inserted into the measuring capsule of a Ksat instrument (UMS, München, Germany).
The Ksat permeameter allows the determination of saturated hydraulic conductivity using two methods, constant-head and falling-head. Both methods are based on the inversion of Darcy's law and fulfil the DIN 19683-9 and DIN 18130-1 standardized procedures [36,37]. Darcy's law defines Ks as a proportionality factor of the amount of water flow through a defined area and the hydraulic gradient.
Ksat allows automated calculation of Ks in the range of 10,000 cm/day down to 0.1 cm/day. In addition, it performs an integrated calculation of Ks at the defined reference temperature according to the dependence of water viscosity on temperature and ensures that there are no water losses due to evaporation during the whole data gathering process.

Aquifer Vulnerability Index and Longitudinal Conductance
The Aquifer Vulnerability Index (AVI) method was developed in Canada by the authors of [38] and uses two variables to quantify a vulnerability index: the thickness of each sedimentary layer above the uppermost saturated aquifer (h) and the estimated hydraulic conductivity of each of these layers (k). The vulnerability index is the sum of the hydraulic resistance (c) of each layer and can be calculated as Equation (1): The k-values for sandy sediments (10 −5 to 10 −1 m/s) are some orders of magnitude higher than those for fine particle size layers (10 −8 to 10 −6 m/s); therefore, hydraulic resistance as defined above is dominated by clayey layers. Hydraulic resistance has the dimension of time (e.g., years) and represents the flux-time per unit gradient for water flowing downward through the layers existing between the surface and the aquifer. The lower the global hydraulic resistance (c), the greater the vulnerability of the underlying aquifer, in absence of preferential flow paths.
This parameter/Equation (1) has the same form as the longitudinal electrical conductance defined by [39] as the second Dar Zarrouk parameter. The Dar Zarrouk parameters were defined to resolve the ambiguity given by the equivalence principle inherent in electrical resistivity interpretation of horizontally layered models, as the parameter is independent of the model chosen. These are easy to compute, and they are related to different combinations of the thickness and resistivity of each geoelectrical layer in the model [16]. For a sequence of n horizontal, homogeneous and isotropic layers of electrical resistivity ρ i and thickness h i the longitudinal conductance, is defined as Equation (2): The relationship between soil parameters (such as clay content, ionic exchange capacity, and vertical hydraulic conductivity) and electrical resistivity enables a vulnerability assessment based on geoelectrical measurements. The results of the measurements can be used to estimate the vertical hydraulic conductivity of the unsaturated zone [40,41]. Generally, clay or fine grain size layers correspond to low resistivities and low hydraulic conductivities, and vice versa. Hence, the protective capacity of the overburden could be considered as being proportional to the ratio of thickness to resistivity-longitudinal conductance (S) [42].
In the present paper, we have calculated the longitudinal conductance parameter from resistivity cross-section data to estimate the protective capacity of the underlying aquifers from percolating contaminants.

Electrical Resistivity Tomography
The 17 ERT cross-sections show resistivity data ranging from 20 Ω·m to more than 3000 Ω·m. The results of the mathematical inversion process have been satisfactory, as the convergence criterion used (root mean square or RMS), has values lower than 4%. From geoelectrical records, three layers can be distinguished according to their resistivity values. The shallowest layer is characterized by values from 80-600 Ω·m and can be identified at the top of cross-sections. The layer has a thickness always identified under 7 m and is interpreted as weathered pyroclasts and areas with transported soils where have been placed ( Figure 2).
Beneath, the geolectrical cross-sections show a layer of fluctuating thickness from 2 to 12 m thickness characterized by values higher than 600 Ω·m. These values are interpreted as porous pyroclasts responses. At this unit, there were significant lateral variations in the resistivity values. The variations reflect a decrease in pyroclast thickness as the distance to the emission center (Pico and Caldera de Bandama) increases.
The third layer is characterized by low resistivity response mainly in the 20-80 Ω·m range and is interpreted as volcanic breccias of the Roque Nublo Group (ignimbritc substrate). Similar outcomes have been obtained by other authors in Tenerife island [43].
The variations in thickness and properties of these three characteristic electrofacies can be clearly seen from the comparison between the cross-sections P1 (Figure 3b) and cross-sections P3 (Figure 3b). The cross-sections were acquired respectively from west to center of the site in the direction of the Bandama Caldera, unveiling an increase of the thickness of the pyroclasts as well as the depth at which the Roque Nublo volcanic debris layer is located.

Hydraulic Conductivity
Hydraulic conductivity results have ranged from minimum values lower than 500 cm/day (consolidated flow pyroclasts) to maximum values above the instrument's measurement limit (20,000 cm/day) for coarse-grained fall pyroclasts (bombs to lapilli). On the other hand, the values of saturated hydraulic conductivity measured by two different

Hydraulic Conductivity
Hydraulic conductivity results have ranged from minimum values lower than 500 cm/day (consolidated flow pyroclasts) to maximum values above the instrument's measurement limit (20,000 cm/day) for coarse-grained fall pyroclasts (bombs to lapilli). On the other hand, the values of saturated hydraulic conductivity measured by two different methods, constant-head and falling-head, have been very congruent, although not equal, as shown in Figure 4a. The interpreted pyroclast layer is narrow and the depth of the low-resistivity ignimbritic substrate (<80 ohm-m) appears at about 5 m deep; (b) ERT cross-section P3. The inferred pyroclast layer is thicker and the low-resistivity ignimbritic substrate is located more than 15 m deep.

Hydraulic Conductivity
Hydraulic conductivity results have ranged from minimum values lower than 500 cm/day (consolidated flow pyroclasts) to maximum values above the instrument's measurement limit (20,000 cm/day) for coarse-grained fall pyroclasts (bombs to lapilli). On the other hand, the values of saturated hydraulic conductivity measured by two different methods, constant-head and falling-head, have been very congruent, although not equal, as shown in Figure 4a.  On average, the values obtained by the constant-head method are between 8 to 30% higher than those obtained by the falling-head method. The authors of [44] consider the constant head method more accurate in the range of hydraulic conductivity between 0.1 to 10 −5 m/s, while the falling-head is better for soils with hydraulic conductivity in the range from greater than 1 m/s to 10 −3 m/s. Since the Bandama Caldera samples cover both groups, it was considered more representative to assign to each sample the arithmetic mean of the two values obtained by each of the methods.
The values of saturated hydraulic conductivity have been compared with the electrical resistivity measured in the golf course itself from the electrical tomography profiles, On average, the values obtained by the constant-head method are between 8 to 30% higher than those obtained by the falling-head method. The authors of [44] consider the constant head method more accurate in the range of hydraulic conductivity between 0.1 to 10 −5 m/s, while the falling-head is better for soils with hydraulic conductivity in the range from greater than 1 m/s to 10 −3 m/s. Since the Bandama Caldera samples cover both groups, it was considered more representative to assign to each sample the arithmetic mean of the two values obtained by each of the methods.
The values of saturated hydraulic conductivity have been compared with the electrical resistivity measured in the golf course itself from the electrical tomography profiles, or by parametric soundings using a Wenner array [45] with 0.2 m of electrode spacing on the same outcrops (Figure 4b).

Longitudinal Conductance
The 18,000 electrical resistivity values from inverted ERT cross-sections were used to estimate the longitudinal conductance (S) value from Equation (2). We have considered h = 20 m-maximum ERT survey penetration depth-and the average of rho values located at the same X and Y position. The Minimum Curvature interpolator was utilized to generate a smooth surface and attempting to honor S data [46]. The Dar Zarrouk parameter S varies from 0.005 Siemens to 5 Siemens. The spatial variation map further infers low S values (0.005-0.02 Siemens) irregularly distributed at the north-eastern, central and southern parts ( Figure 5). S values greater than 0.1 Siemens were mainly identified in the central and southern sectors. The results show the highest resolution in areas with ERT data.
The protective capacity is assumed to be directly proportional to the longitudinal conductance (S). Accordingly, the overburden protective capacity was evaluated using the total longitudinal unit conductance (S). In the studied area lower S values generally indicate a relatively weak succession of fine grain size sediments overburden together with greater proximity to the emission center of the eruption and are given the highest priority in terms of aquifer protection studies as it implies the potential infiltration of contaminants into the aquifer [47]. at the same X and Y position. The Minimum Curvature interpolator was utilized to generate a smooth surface and attempting to honor S data [46]. The Dar Zarrouk parameter S varies from 0.005 Siemens to 5 Siemens. The spatial variation map further infers low S values (0.005-0.02 Siemens) irregularly distributed at the north-eastern, central and southern parts ( Figure 5). S values greater than 0.1 Siemens were mainly identified in the central and southern sectors. The results show the highest resolution in areas with ERT data. The protective capacity is assumed to be directly proportional to the longitudinal conductance (S). Accordingly, the overburden protective capacity was evaluated using the total longitudinal unit conductance (S). In the studied area lower S values generally indicate a relatively weak succession of fine grain size sediments overburden together with greater proximity to the emission center of the eruption and are given the highest priority in terms of aquifer protection studies as it implies the potential infiltration of contaminants into the aquifer [47].

Discussion
The hydraulic conductivity of volcanic formations is a difficult parameter to measure and usually presents a high anisotropic ratio causing the infiltrating water to prefer the horizontal flow component while the vertical flow remains as a secondary path. Moreover, it has a wide variability due to genetics, petrochemical composition and geological history, including deposition mechanisms, alteration, lithification or the existence of subsequent fractures and compactions. Consistently, the hydraulic conductivity of volcanic formations is expressed in wider ranges of values than in other formations [48]. Table 1 presents hydraulic conductivity of main volcanic formations of the study area obtained by usual hydrodynamic technics. In general, young and non-welded pyroclasts have high permeability and altered or consolidated pyroclasts have low values [49]. Hydraulic conductivity values by Ksat equipment are consistent with these cited wide range (Table 1). Table 1. Hydraulic conductivity values obtained by previous studies [49], and maximum and minimum K values measured by Ksat equipment for this study in volcanic formations of the zone. Hydraulic conductivity could be estimated indirectly from electrical resistivity values [50]. Nevertheless, this correlation must be made based on local tests and with reservations, since electrical resistivity is also a function of the degree of saturation and the electrical conductivity of the soil water. On golf courses, if resistivity measurements are made after the irrigation procedure with an excess of water, the subsoil can be considered to have a moisture content close to field capacity. According to the authors of [51], the electrical resistivity values tend asymptotically to the saturation value under these conditions both in pyroclastics volcanic soils (Figure 6a) and in volcanic soils with ignimbrites ( Figure 6b). Table 1. Hydraulic conductivity values obtained by previous studies [49], and maximum and minimum K values measured by Ksat equipment for this study in volcanic formations of the zone.

Volcanic Formation Horizontal Hydraulic Conductivity (m/d) K (m/d) Obtained in This Study
Recent basalts 5-1000 200 for coarse-grained pyroclasts fall deposits Volcanic breccias (Roque Nublo Group) <0.002-0.5 <5 m/d for consolidate pyroclasts flow deposits Hydraulic conductivity could be estimated indirectly from electrical resistivity values [50]. Nevertheless, this correlation must be made based on local tests and with reservations, since electrical resistivity is also a function of the degree of saturation and the electrical conductivity of the soil water. On golf courses, if resistivity measurements are made after the irrigation procedure with an excess of water, the subsoil can be considered to have a moisture content close to field capacity. According to the authors of [51], the electrical resistivity values tend asymptotically to the saturation value under these conditions both in pyroclastics volcanic soils (Figure 6a) and in volcanic soils with ignimbrites (Figure 6b). The preferential infiltration zones have been delimited by correlating electrical resistivity and vertical hydraulic permeability of the different geological units of the unsaturated zone that outcrop at the east edge of the golf course (represented as recent volcanoes and Roque Nublo Group in Figure 1b). The pyroclastic layers with the highest electrical resistivity have the highest porosity and, in turn, the ones with the highest hydraulic permeability. Transported soil was measured in P5, P6, P8, P14, P15 and P16 (Figure 7). Those soils present variable hydraulic properties due to their structure and content in organic The preferential infiltration zones have been delimited by correlating electrical resistivity and vertical hydraulic permeability of the different geological units of the unsaturated zone that outcrop at the east edge of the golf course (represented as recent volcanoes and Roque Nublo Group in Figure 1b). The pyroclastic layers with the highest electrical resistivity have the highest porosity and, in turn, the ones with the highest hydraulic permeability. Transported soil was measured in P5, P6, P8, P14, P15 and P16 (Figure 7). Those soils present variable hydraulic properties due to their structure and content in organic matter and will probably be less deep which explains the need to add transported soil [19]. Moreover, in this study, its narrow thickness seems not to modify the global average of electrical resistivity values as we use 20 m for the resistivity assessment presented in Figure 7. Conversely, we could identify a reduction of resistivity values in the closest part to the emission point (northeast part) where the high resistivity layer of pyroclasts is placed deeper (more than 7 m). The existence of water oozes under the lower and thicker soil layer and a water gallery in the slope of the Caldera in fractured ignimbrite under the pyroclastic layer corroborates the results [6]. [19]. Moreover, in this study, its narrow thickness seems not to modify the global average of electrical resistivity values as we use 20 m for the resistivity assessment presented in Figure 7. Conversely, we could identify a reduction of resistivity values in the closest part to the emission point (northeast part) where the high resistivity layer of pyroclasts is placed deeper (more than 7 m). The existence of water oozes under the lower and thicker soil layer and a water gallery in the slope of the Caldera in fractured ignimbrite under the pyroclastic layer corroborates the results [6]. Figure 7. Distribution of the shallow subsurface electrical resistivity interpolated using kriging algorithm [52] from ERT profiles (dark grey lines) acquired over the golf course. Dotted areas represent the location of transported soils.
The pollution of groundwater as a result of different antropogenic activities, including the irrigation of golf courses with an excess of reclaimed water, is one of the main obstacles faced by most of the administrations of the European Union member states to achieve the objectives of the Water Framework Directives [53]. To this end, it is essential to assess the best reclaimed water irrigation management practices based on the vulnerability to contamination. It is also necessary to take extreme precautions in vulnerable areas. Both concepts are based on a better knowledge of the infiltration and migration of contaminants through the unsaturated zone and the necessity for defining the protective Figure 7. Distribution of the shallow subsurface electrical resistivity interpolated using kriging algorithm [52] from ERT profiles (dark grey lines) acquired over the golf course. Dotted areas represent the location of transported soils.
The pollution of groundwater as a result of different antropogenic activities, including the irrigation of golf courses with an excess of reclaimed water, is one of the main obstacles faced by most of the administrations of the European Union member states to achieve the objectives of the Water Framework Directives [53]. To this end, it is essential to assess the best reclaimed water irrigation management practices based on the vulnerability to contamination. It is also necessary to take extreme precautions in vulnerable areas. Both concepts are based on a better knowledge of the infiltration and migration of contaminants through the unsaturated zone and the necessity for defining the protective properties naturally occurring in geologic layers. The variation in reclaimed water quality through time also supports the use of vulnerability models. As previous studies demonstrated the desalination treatment implemented in 2002, reduced the electrical conductivity of irrigation water from 2800 to 1000 µS/cm, affecting the infiltration soil rate of the golf course [19].
Of particular significance is defining vertical travel times (TTs) through layers located above aquifers to prevent contamination from cultural activities. Surface geoelectric provide a fast and economical field method that can be used to assess the protective properties of geologic layers. In particular, the TT through unsaturated layers is theoretically linearly related to the longitudinal unit conductance (S) of the layers with an estimated standard deviation of 2.9 years by authors of [54].
Nevertheless, it must be considered that the longitudinal conductivity model is considered a semiquantitative assessment and requires a site-specific classification to rate the protective capacity of the unsaturated area [55,56].
We have been followed the criteria of the AVI methodology to assign the vulnerability categories. TTs of more than 3 years have been identified in the southern and eastern zones of the studied site, these being, a priori, the areas most protected from surface contamination of the studied area ( Figure 8). On the other hand, transit times of less than 1 year, and therefore areas vulnerable to surface contamination, are located mainly in the north and eastern parts.
properties naturally occurring in geologic layers. The variation in reclaimed water quality through time also supports the use of vulnerability models. As previous studies demonstrated the desalination treatment implemented in 2002, reduced the electrical conductivity of irrigation water from 2800 to 1000 µ S/cm, affecting the infiltration soil rate of the golf course [19].
Of particular significance is defining vertical travel times (TTs) through layers located above aquifers to prevent contamination from cultural activities. Surface geoelectric provide a fast and economical field method that can be used to assess the protective properties of geologic layers. In particular, the TT through unsaturated layers is theoretically linearly related to the longitudinal unit conductance (S) of the layers with an estimated standard deviation of 2.9 years by authors of [54].
Nevertheless, it must be considered that the longitudinal conductivity model is considered a semiquantitative assessment and requires a site-specific classification to rate the protective capacity of the unsaturated area [55,56].
We have been followed the criteria of the AVI methodology to assign the vulnerability categories. TTs of more than 3 years have been identified in the southern and eastern zones of the studied site, these being, a priori, the areas most protected from surface contamination of the studied area ( Figure 8). On the other hand, transit times of less than 1 year, and therefore areas vulnerable to surface contamination, are located mainly in the north and eastern parts. In the case of the study area, the data provided in this work, conveniently crosschecked with geological data from the Caldera wall, will allow its calculation to be included in future models of contaminant transport through the unsaturated zone. In the case of the study area, the data provided in this work, conveniently crosschecked with geological data from the Caldera wall, will allow its calculation to be included in future models of contaminant transport through the unsaturated zone.

Conclusions
The research results show that electrical resistivity tomography is a suitable technique to investigate quickly and non-destructively the geometry and lithological characteristics of the subsoil and to assess the best reclaimed water irrigation management practices and the vulnerability to contamination of groundwater beneath golf courses, even in complex geological environments, as in the case of the Bandama Golf Course.
The electrical resistivity values have made it possible to identify each of the lithological units that make up the subsoil of the golf course, providing a general model that agrees with the edaphological observations made based on the geological knowledge of the volcanic structure where the course is located. The model obtained provides detailed information on the lateral and vertical variability of each of the layers and, based on an empirical correlation between the values of electrical resistivity and hydraulic permeability, makes it possible to delimit the preferential zones of subsurface drainage that may represent a greater risk to the vulnerability of the underlying aquifer.
The AVI method is a quantitative method that allows determining vulnerability in terms of the transit time of the contaminant through the unsaturated zone. The transit time can be estimated by indirect methods, based on the information provided by electrical resistivity tomography without affecting the game development and preserving the playground integrity.