Geophysical Characterization of Hydraulic Properties around a Managed Aquifer Recharge System over the Llobregat River Alluvial Aquifer (Barcelona Metropolitan Area)

Managed aquifer recharge using surface or regenerated water plays an important role in the Barcelona Metropolitan Area in increasing storage volume to help operators cope with the runoff variability and unexpected changes in surface water quality that are aggravated by climate change. The specific aim of the research was to develop a non-invasive methodology to improve the planning and design of surface-type artificial recharge infrastructures. To this end, we propose an approach combining direct and indirect exploration techniques such as electrical resistivity tomography (ERT), frequency domain electromagnetics and data from double-ring infiltration tests, trial pits, research boreholes and piezometers. The ERT method has provided much more complete and representative information in a zone where the recharge project works below design infiltration rates. The geometry of the hydrogeological units and the aquifer-aquiclude contact are accurately defined through the models derived from the interpretation of ERT cross-sections in the alluvial aquifer setting. Consequently, prior to the construction of recharge basins, it is highly recommended to conduct the proposed approach in order to identify the highest permeability areas, which are, therefore, the most suitable for aquifer artificial recharge.


Introduction
The Barcelona Metropolitan Area has a drinking water demand for almost three million inhabitants, although its sources of water are scarce because of its Mediterranean climate. Their main local sources are surface water from the Llobregat River and groundwater from the underlying aquifer. Many European rivers have fairly constant flows > 1500 m 3 /s, but Mediterranean rivers are characterized by their low average flows with intermittently high flow peaks, such as the Llobregat River (22 m 3 /s of annual average flow but with a 2-5 m 3 /s flow most days of the year). In addition to this water scarcity, the management of drinking water in the Mediterranean area has, consequently, some other extra difficulties in comparison with other areas in Europe, such as the industrial impacts, the pollution, the wastewater management, changes in land use and the demand of a minimum ecological flow in the context of the Water Framework Directive. For these reasons, good management integrating all possible sources of water in the Barcelona Metropolitan Area is required [1].
Among other good water management practices, the artificial recharge of aquifers is perhaps the most effective in regulating water resources and ensuring water supply during dry seasons [2]. Managed aquifer recharge (MAR) is the intentional recharging of water to suitable aquifers for subsequent recovery or to achieve environmental benefits: the managed process ensures adequate protection of human health and the environment [3]. MAR is a widely employed technique [4] and using both surplus water and regenerated water can be a good alternative in densely populated areas with high water demand, where the extraction of this resource might exceed the natural inflow to the aquifer [5]. In the case of the Barcelona Metropolitan Area, MAR plays an important role in increasing storage volume to help operators cope with the runoff variability in Barcelona catchments and unexpected changes in surface water quality that are aggravated by climate change. Historically, public administration, the Catalan Water Agency (ACA) and water supply operators (the SGAB Company) have applied several techniques to artificially improve the groundwater resources-using surface water or regenerated water-in the Low Llobregat Valley and in the Delta of the Llobregat River, such as injection wells, river bed scarification or infiltration ponds [6].
Regardless of the purpose of their use, accurate knowledge of the geology and hydrogeology of the medium is strongly required, both in the planning and operational stages, in order to better understand the hydrological response of the MAR system [7]. Particularly, infiltration ponds' systems usually require permeable surface soils to get high infiltration rates and to minimize land use requirements. The non-saturated zone (NSZ) should be free of fine/textured materials such as clays that unduly restrict downward flow and form perched groundwater that waterlogs the recharge area and reduces infiltration rates [3]. Ubiquitously, insufficient site characterization and soil heterogeneities render prediction of an infiltration pond's performance fundamentally uncertain and are a determining factor for optimal decision making in maintenance operations [8,9].
The use of remote sensing techniques such as satellite images together with Geographic Information Systems (GIS) was successfully implemented to evaluate favorable recharge areas [10][11][12]. However, they are usually constrained as a preliminary reference in selecting suitable sites for groundwater resources management and in narrowing down the target areas for conducting detailed hydrogeological surveys at the catchment scale [13].
Classical detailed hydrological approaches such as quantifying hydraulic conductivity from pumping and recovery tests require the existence of available wells and piezometers. In addition, such tests provide only punctual assessments of the aquifer parameters and hydraulic properties. Therefore, considering the amount of data required to adequately characterize a NSZ around a MAR project at the necessary scale can be cost-prohibitive and functionally impractical due to natural subsurface heterogeneity [14].
The potential benefits of including geophysical data in hydrogeological site characterization have been stated in numerous studies [15][16][17][18][19]. Among others, geophysical measurements are often spatially exhaustive and can be acquired at relatively low cost, making them very useful for supplementing conventional measurements. Electrical resistivity tomography (ERT) and frequency domain electromagnetics (FDEM) are geophysical methods particularly useful for studies of hydrologic processes in the near-surface because of the high correlation between electrical conductivity and water saturation in the vadose zone [20][21][22][23][24][25][26][27], for surface MAR facilities monitoring [28][29][30] and for geological heterogeneities characterization in alluvial aquifers [31]. When geophysical measurements are used in an early phase of site investigations, they may provide useful input to the placement of new research boreholes or wells. This would potentially limit the number of wells, and information from drilling operations can be more valuable [32].
The specific aim of our research was to develop a non-invasive methodology to improve the planning and design of surface-type artificial recharge infrastructures. To this end, we propose an approach combining direct and indirect exploration techniques in four phases. The first two phases, which correspond to the geophysical characterization of the NSZ and the saturated zone (SZ), provides a detailed characterization of the subsurface of the Ca n'Albareda infiltration pond in terms of the electrical resistivity parameter, using ERT and FDEM and data from trial pits, boreholes and piezometers. The last two phases (infiltration rate and assessment of artificial recharge potential) have been applied to the Ca n'Albareda alluvial aquifer. ERT and infiltration tests were carried out to obtain a local relationship between the measured physical property (electrical resistivity) and the one of interest for the project (infiltration rate) and highlight the lateral and vertical heterogeneities of the subsurface and the existence of clayey sediments.

Site Description
The study areas are the artificial recharge ponds and the meander of Ca n'Albareda. They are located northeast of the Iberian Peninsula within the Metropolitan Area of the city of Barcelona, on the alluvial plain of the Llobregat River and about 25 km upstream from its mouth to the Mediterranean Sea ( Figure 1). The Llobregat River is 168.5 km long and is considered one of the most studied and monitored river systems in Europe [33]. However, direct and quantitative hydrogeological data from the study zone were virtually non-existent and unrepresentative of the expected heterogeneity of the geology and hydraulic conductivity in a sedimentary environment such as Ca n'Albareda meander [34].
Water 2020, 12, x FOR PEER REVIEW 3 of 22 of the electrical resistivity parameter, using ERT and FDEM and data from trial pits, boreholes and piezometers. The last two phases (infiltration rate and assessment of artificial recharge potential) have been applied to the Ca n'Albareda alluvial aquifer. ERT and infiltration tests were carried out to obtain a local relationship between the measured physical property (electrical resistivity) and the one of interest for the project (infiltration rate) and highlight the lateral and vertical heterogeneities of the subsurface and the existence of clayey sediments.

Site Description
The study areas are the artificial recharge ponds and the meander of Ca n'Albareda. They are located northeast of the Iberian Peninsula within the Metropolitan Area of the city of Barcelona, on the alluvial plain of the Llobregat River and about 25 km upstream from its mouth to the Mediterranean Sea ( Figure 1). The Llobregat River is 168.5 km long and is considered one of the most studied and monitored river systems in Europe [33]. However, direct and quantitative hydrogeological data from the study zone were virtually non-existent and unrepresentative of the expected heterogeneity of the geology and hydraulic conductivity in a sedimentary environment such as Ca n'Albareda meander [34].  [35]. Electrical resistivity tomography (ERT) locations presented in this paper are quoted in light red. The meander covers an area of about 60 hectares and comprises part of the smallest and largest bed of the river (35 masl) and its right (42 masl) and left (39 masl) banks. The ponds of Ca n'Albareda are excavated structures within the meander, have a longitudinal orientation east-west and are located on the right bank of the Llobregat River. They form a new wetland area composed of two ponds, each with a differentiated function. The first is the decanting pond, located upstream, where the suspended matter from the water is retained prior to passage to the second, and the second is the infiltration pond, where the water recharge process takes place. The infiltration pond has an average depth of four meters, a total length of 220 m and an infiltrating surface of 6000 m 2 . The operating procedure chosen by facility managers was to maintain the water level within the pond from 0.6 to 1.2 m above the ground. It was estimated at a 1-2 m 3 /m 2 ·day infiltration rate and a total volume of about 1.6 × 10 6 m 3 /year, considering 180 annual working days.
The Ca n'Albareda sites are located in the hydrogeological unit of the Cubeta de Sant Andreu de la Barca basin, where the Llobregat River runs following a directional fault path. The fault generates a great asymmetry between the two banks of the river, separating Paleozoic schists and Mesozoic limestones to the west from the east Cainozoic sediments.
The geological framework of the unit was generated by the excavation and sedimentation of the river itself, producing the formation of different river terraces, ranging from the oldest, but located at an upper-level terrace (Qt 3 ) to the current river path (Qt 0 ). River terraces are composed of gravel, sand and silt sediments and constitute a phreatic alluvial aquifer of tabular geometry. The maximum thickness is about 30 m and the outcrop extension ranges between eight and ten km 2 [35]. The hydraulic conductivity is higher than 100 m/day, the storage coefficient ranges from 0.1 to 0.3 and it has estimated groundwater withdrawn of 6 × 10 6 m 3 /year [36].
Clays and silts deposited during Miocene are located nine meters below the infiltrating surface ( Figure 2). They are described as red claystone, present hydraulic conductivity values from 10 −5 to 10 −3 m/day and act as an aquiclude [35]. The area is characterized by a Mediterranean climate with a warm thermal regime in summer and moderate cold in winter. The annual rainfall is on average between 600 and 650 mm/year. Maximum rain values are recorded in autumn and the driest month is July, coinciding with the highest potential evapotranspiration values [37].  The area is characterized by a Mediterranean climate with a warm thermal regime in summer and moderate cold in winter. The annual rainfall is on average between 600 and 650 mm/year. Maximum rain values are recorded in autumn and the driest month is July, coinciding with the highest potential evapotranspiration values [37].

Geophysical Surveys
Two complementary geophysical methods were applied in order to perform the characterization of lateral and vertical lithological variations of the SZ below the infiltration pond: ERT and FDEM.
ERT was used to characterize NSZ, to obtain a relationship between electrical resistivity and infiltration rate and to characterize the recharge potential in the meander of Ca n'Albareda among seven geophysical acquisition campaigns carried out in the 2010-2013 period ( Table 1). ERT data was acquired with a Syscal Pro device (IRIS instruments, Orléans, France). The system features an internal switching board for 48 electrodes and an internal 250 W power source. A Wenner-Schulmberger array was chosen because it is moderately sensitive to both horizontal and vertical structures and has a relatively good signal strength [38].
The physical parameter measured by the ERT method is the apparent electrical resistivity of the subsurface. The software used for the inversion of the two-dimensional (2D) ERT data was RES2DInv v 3.54, [39] and RES3DInv 2.14 [40] and ERTlab [41] were used for the inversion of three-dimensional (3D) ERT data. The subsurface is divided into cells of fixed dimensions and the inversion procedure is based on the smoothness-constrained least-squares method. The resistivities are adjusted iteratively until a satisfactory agreement between the input data and the model responses is achieved, based on a nonlinear optimization technique by least-squares fitting [42]. During the inversion process, the root-mean-square value of the difference between experimental data and the updated model response is used as a criterion to assess the convergence.
In the present paper, the robust least-squares method was selected, after making a comparison with the smoothing constraint method in SZ site characterization. The method assumes that the subsurface consists of a few homogeneous regions with a sharp interface between them. Such an inversion scheme is the logical choice where the subsurface comprises units with sharp boundaries in order to accurately determine both layer boundary locations and layer resistivities. Indeed, it produces models by minimizing the absolute value of data misfit, making it more efficient in removing noise compared to other inversion methods [43].
FDEM apparent electrical conductivity data were measured using a EM31-MK2 instrument (Geonics Limited, Ontario, Canada). The system consists of transmitting and receiving coils spaced 3.66 m apart and operates at a frequency of 9.8 kHz. When operated at a low induction number, i.e., where the electromagnetic skin depth greatly exceeds the transmitter/receiver coil space, the instrument provides the apparent, or depth-averaged, conductivity at each measurement position [44]. In our case, two ground conductivity measurements were performed on each station to validate the ERT results on the SZ characterization survey. The coils can be positioned either in a vertical (VD) or horizontal (HD) dipolar configuration. The horizontal configuration was used to investigate the first three meters and the vertical configuration the first six meters.

Infiltration Rate
The effectiveness of artificial recharge ponds is controlled in large part by the infiltration capacity of the host topsoil [8]. The infiltration capacity of soil is the maximum amount of water that can absorb in a unit of time and with defined soil conditions, being able to define as the rate at which water penetrates inside the soil through its surface. The infiltration rate is equivalent to the vertical component of hydraulic conductivity and basically depends on the texture and structure of the soil and can be determined from in-situ measurements or runoff analysis in small basins [45].
To obtain the infiltration rate values, nine in-situ infiltration tests were performed with a double-ring or Munz infiltrometer [46] at the midpoint of the parametric 2D ERT cross-section (23.5 m long). The double-ring infiltrometer is an instrument that provides direct information on the subsurface infiltration rate, or field-saturated hydraulic conductivity, in a simple way [47]. Nevertheless, the test needs a water supply, is demanding in terms of time spent and is usually only representative of the tested point and the shallowest part of the subsurface.
To carry out the infiltration tests, different points representative of the sediments' outcropping in the Cubeta were chosen ( Table 2). In the case of the Ca n'Albareda recharge pond, a few months after reopening in 2010, the infiltration ponds operated at significantly lower infiltration rates than design parameters (1-2 m 3 /m 2 ·day). A pilot excavation was placed in a pond sector to improve the system and two infiltration tests were carried out using the same procedure: one in the original surface (It 1 ) and the other over the excavated area (It 2 ). The two tests' results have already been used in previous research on monitoring of clogging process and maintenance actions outcomes [30].

Assessment of Aquifer Recharge
Aquifer transmissivity is an important issue in groundwater management projects because it illustrates the rate of groundwater movement throughout a unit area in the aquifer. By definition, it is the "ability of the aquifer to transmit groundwater throughout its saturated thickness" [48] and it is calculated by multiplying hydraulic conductivity and the saturated thickness in a homogeneous media [45]. Areas with high aquifer transmissivity values have been considered in the literature favorable for artificial groundwater recharge [49,50].
Hydraulic transmissivity values are generally determined from pumping tests. Often, hydraulic transmissivity values are average values in a large aquifer volume, can have large variations from one point to another of the aquifer, and in many cases, pumping tests are not easy to perform and are generally relatively costly in time and resources. Therefore, while always necessary for good hydrogeological characterization, it is very difficult to obtain a representative distribution of transmissivity values of an entire aquifer unit.
In this study, we used the information obtained with geoelectric methods as they have a much more representative spatial distribution of the area in order to identify both the most and least favorable areas for the exploitation of groundwater resources and to perform artificial recharge operations in the Ca n'Albareda meander. However, geoelectrical data inversion and interpretation are unavoidably limited by the ambiguity expressed in the principle of equivalence, which states that an experimental resistivity curve can be explained as having been produced by many physically equivalent models [51]. This ambiguity mainly affects the individual assessment of the parameters of each layer, but not the geoelectric properties of the model assembly. In particular, the use of the combination of layer resistivity and thickness (h) in the so-called Dar Zarrouk parameters S and T may be of direct use for the evaluation of hydrologic properties of aquifers [52].
The T Dar Zarrouk parameter or transverse resistance-defined by the product of h and electrical resistivity of a layer-can be correlated with their hydraulic transmissivity [53][54][55]. Negative statistical relationships have been reported in aquifers where water circulates preferentially by fractures (karst type) [56], and positive in granular aquifers, as in our case [57]. In this paper, we have calculated the T Dar Zarrouk parameter, from resistivity cross-section data, using a similar scheme as aquifer transmissivity calculation: from the water table to the bottom of the aquifer unit.

Geophysical Characterization of the Saturated Zone (SZ)
The design of geophysical surveys and the selection of the multi-electrode arrays was planned considering the depth of investigation, the length and resolution required and the expected structure derived from hydrogeological background knowledge. A lithological log of a borehole drilled previously between the right bank of the Llobregat River and the infiltration pond of Ca n'Albareda (S2), and the logs of two additional piezometers (CA1 and CA2) equipped inside the infiltration pond, in parallel to the geophysical survey, have been an invaluable support for our research ( Figure 2).
In general, geological logs from boreholes showed a coarsening upward trend. Two main hydrogeological units could be easily identified: an upper unit mainly constituted by gravels and brown sands (aquifer) and a bottom unit, mainly formed by red claystone (low permeability layer or aquiclude). The contact between the two formations was located at 9.5 m below the surface of the infiltration pond.
In order to characterize the geometry of the aquifer below the recharge pond and the possible geological discontinuities of the subsurface, the SZ1 and SZ2 2D ERT cross-sections were acquired ( Figure 1) and a 40 × 30 m 3D ERT array has been built in the center of the infiltration pond. The 2D ERT cross-sections of 235 m allowed to reach a research-depth close to 50 m and a resolution of one point every five meters in both directions.
Two inversion methods have been used for each of the two cross-sections: smoothing-constrained and robust. In all cases, the results of the mathematical inversion process have been satisfactory, as the convergence criterion used (root mean square or RMS) has values always below 5%. However, after comparing the two sections obtained in each case, the robust method was chosen due to more consistent geological interpretation and better definition of the geometry of the identified lithological units, especially the contact between the aquifer and the aquiclude.
The subsequent subsurface characterization using electrical conductivity or resistivity depends on several factors, such as soil water content, grain size distribution, porosity and permeability. For instance, air-filled void soil type will have higher geoelectrical resistivity values contrary to a water-filled void soil type. Also, in fine-grained subsoil materials such as clay, where soil water content is higher, the observed geoelectrical resistivity is usually low [58].
In our case, the electrical resistivity values obtained from the ERT cross-sections ranged between 5 and 200 Ωm. A general trend of increase from the bottom to the top was observed, with a notable variation in upper-level values. Three layers can be distinguished in all cross-sections according to their electrical resistivity values (Figure 3). At the top of the SZ1 ERT cross-section, the electrical resistivity values are relatively high (>200 Ωm), with an increase in electrical resistivity values from position 120 m to the end of the profile. In the intermediate layer, an area of low electrical resistivity values is identified at the bottom of the profile, interpreted as clays. The geometry of the contact between clean gravels and gravels with a clayey matrix is quite regular and located at an average depth of about 11 m. In cross-section SZ2, the geometry of the contact between the layer of gravel and the layer of clay substrate is less regular and located about 10-13 m deep. A variation is also identified at the upper part of the subsurface, with higher resistivity values towards both sides of the profile, relatable to textural variations in the gravel-sized sediments (Figure 3). At the bottom of the cross-section, the distribution of electrical resistivity values is almost homogeneous (20-50 Ωm). A 3D resistivity model (40 × 30 m) of the central part of the recharge pond subsoil has been obtained from 3D ERT acquisition using the ERTLab software ( Figure 4). The electrical resistivity values in the model range from 64 to 662 Ωm. The investigation depth was eight meters deep and allows us to define in detail the variability of the topmost part, where the gravels with the clayey matrix are present. A gradual increase appears, from the base to the top, of the resistivity values, and variations in the upper levels, emphasizing three sectors with resistivity values above 400 Ωm. In cross-section SZ2, the geometry of the contact between the layer of gravel and the layer of clay substrate is less regular and located about 10-13 m deep. A variation is also identified at the upper part of the subsurface, with higher resistivity values towards both sides of the profile, relatable to textural variations in the gravel-sized sediments (Figure 3). At the bottom of the cross-section, the distribution of electrical resistivity values is almost homogeneous (20-50 Ωm). part of the subsurface, with higher resistivity values towards both sides of the profile, relatable to textural variations in the gravel-sized sediments (Figure 3). At the bottom of the cross-section, the distribution of electrical resistivity values is almost homogeneous (20-50 Ωm). A 3D resistivity model (40 × 30 m) of the central part of the recharge pond subsoil has been obtained from 3D ERT acquisition using the ERTLab software ( Figure 4). The electrical resistivity values in the model range from 64 to 662 Ωm. The investigation depth was eight meters deep and allows us to define in detail the variability of the topmost part, where the gravels with the clayey matrix are present. A gradual increase appears, from the base to the top, of the resistivity values, and variations in the upper levels, emphasizing three sectors with resistivity values above 400 Ωm.

Geophysical Characterization of the Non-Saturated Zone (NSZ)
The characterization of the NSZ was performed in two sub-phases or stages and has been validated with the direct information provided by trial pits excavated at different points of the infiltration pond. The first stage was carried out using both ERT and FDEM devices, whereas for the second stage, only the ERT equipment was used.
In the first stage, a point-to-point survey inside the infiltration pond was designed in the form of four rectangular grids, with measurements separated five meters from each other. Then, the apparent electrical conductivity values obtained with the EM-31 instrument have been interpolated using the Radial Basis Function method [59], with a spacing of five meters and limiting the interpolation zone to the edges of the infiltration pond. To compare with the results of the EM-31 survey, a 235 m long 2D ERT was acquired on the south margin of the recharge pond ( Figure 5).
Due to FDEM's higher acquisition speed per surface unit compared to ERT, a qualitative conductivity map has been built first). The electrical conductivity map-obtained from FDEM data-shows zones with low conductivity (7-10 mS/m) located at the western part of the pond. These values are increasing within the east direction and maximum values were identified at the eastern part of the pond (13-15 mS/m). The results of the ERT cross-section corroborate this trend as the area of the pond with lower electrical conductivity values, the 2D cross-section, shows the highest resistivity values at the same depth. Both results were used to plan the second-and more detailed-characterization of the NSZ using the highest resolution technique: ERT. It is in that central and western area with lower electrical conductivity values in which a more detailed survey using 23.5 m long 2D ERT was performed.
In the second stage, 11 ERT cross-sections of 47 m in length (9 m of maximum depth) and 38 ERT cross-sections of 23.5 m in length (4.5 m deep) were obtained. Almost 20,000 apparent resistivity values from the 38 high-detailed ERT acquisitions were used to build a 74 × 23.5 × 6.09 m 3D resistivity model. RES3DINV software allows to carry a 3D inversion-in that the resistivity values can vary in all three directions simultaneously during the inversion process-and plot results in resistivity depth-slices.
Both 2D and 3D models show heterogeneous resistivity values below the surface of the infiltration pond. Three layers are easily identified according to their different electrical resistivity values: On the other hand, the Association of Water Users of the Sant Andreu de la Barca basin (CUACSA) carried out a characterization of the basement of the infiltration pond. Eight mechanical trenches were excavated, with an offset of approximately 25 m from each other. The result of the trial pits shows the existence of two distinct geological layers: an upper layer consisting mainly of a gravel-sized material with a silty-clay matrix and a lower unit consisting mainly of coarse gravel and sand. The thickness of the upper unit ranges from 0.6 (pit N4) to 1.3 m (pit N6) ( Figure 6). The thickness of the lower unit cannot be defined from the mechanical trenches because its low boundary is placed below the length of the arm of the excavator used. Due to FDEM's higher acquisition speed per surface unit compared to ERT, a qualitative conductivity map has been built first). The electrical conductivity map-obtained from FDEM datashows zones with low conductivity (7-10 mS/m) located at the western part of the pond. These values are increasing within the east direction and maximum values were identified at the eastern part of the pond (13-15 mS/m). The results of the ERT cross-section corroborate this trend as the area of the pond with lower electrical conductivity values, the 2D cross-section, shows the highest resistivity In this context, the ERT cross-sections might show the morphology and thickness variations of the alluvial gravel unit and the shallowest layer characterized by its clayey matrix. The three geoelectric layers described above can therefore be interpreted from top to bottom, as clay matrix gravel, coarse gravels with sand matrix and water-saturated gravels (Figure 7). It should be considered that the base layer has been interpreted using the water table position measured directly at the piezometer CA2 during ERT acquisitions. The water table is located at a depth of 3.4 below the surface of the pond.


A discontinuous layer with low values located between 0 and 1.5 m deep  A layer with relatively high resistivity between 0.6 and 3 m deep  A layer with lower resistivity values located at the bottom of all the cross-sections On the other hand, the Association of Water Users of the Sant Andreu de la Barca basin (CUACSA) carried out a characterization of the basement of the infiltration pond. Eight mechanical trenches were excavated, with an offset of approximately 25 m from each other. The result of the trial pits shows the existence of two distinct geological layers: an upper layer consisting mainly of a gravelsized material with a silty-clay matrix and a lower unit consisting mainly of coarse gravel and sand. The thickness of the upper unit ranges from 0.6 (pit N4) to 1.3 m (pit N6) ( Figure 6). The thickness of the lower unit cannot be defined from the mechanical trenches because its low boundary is placed below the length of the arm of the excavator used. In this context, the ERT cross-sections might show the morphology and thickness variations of the alluvial gravel unit and the shallowest layer characterized by its clayey matrix. The three geoelectric layers described above can therefore be interpreted from top to bottom, as clay matrix gravel, coarse gravels with sand matrix and water-saturated gravels (Figure 7). It should be considered that the base layer has been interpreted using the water table position measured directly at the piezometer CA2 during ERT acquisitions. The water table is located at a depth of 3.4 below the surface of the pond.

Infiltration Rate
An empirical relationship has been derived from the electrical resistivity values and infiltration rates obtained at all the nine test sites (Figure 8). The relationship aims to obtain an estimate of the infiltration rate (Ir) based on the central and shallowest data of the parametric ERT. The duration of each of the tests ranged from two hours from point A3 to five hours from point A7, and the amount

Infiltration Rate
An empirical relationship has been derived from the electrical resistivity values and infiltration rates obtained at all the nine test sites (Figure 8). The relationship aims to obtain an estimate of the infiltration rate (Ir) based on the central and shallowest data of the parametric ERT. The duration of each of the tests ranged from two hours from point A3 to five hours from point A7, and the amount of water to obtain the required saturation conditions ranged from the 60 L of the point A7 to 300 L of point A3.
Water 2020, 12, x FOR PEER REVIEW 14 of 22 0.9, but the positive slope of the relationship suggests that areas with ERT high values will have a high vertical hydraulic conductivity (k).

Assessment of Recharge Potential
The results of the 23 ERT cross-sections used for the characterization of the subsoil below the Ca n'Albareda meander ( Figure 1) show a high contrast in the electrical resistivity values, ranging from 5 to 1300 Ωm. This fact demonstrates the variability of the lithology and the hydraulic properties of sediments existing in the subsurface. Nevertheless, there is an overall trend of an increase in the electrical resistivity values from the bottom to the top and a noticeable variation in the upper levels.
By examining each ERT cross-section one by one, a layer with relatively low electrical resistivity values can be identified between the surface and about seven meters depth, and an intermediate layer with higher values. At the bottom, a low resistivity layer is identified, and it can be considered the impervious basement of the alluvial aquifer. The depth of this limit in the ERT cross-sections ranges between 15 (Rp1) and almost 40 m (Rp7) and, with the meander being a topographically flat area, defines an aquifer with a very variable thickness depending on the part of the meander considered. In addition, as can also be seen at the bottom of Figure 9, the electrical resistivity response of the materials forming the aquifer is different due to the heterogeneity of alluvial sediments. The results of the infiltration tests show higher Ir in higher resistivity areas and lower Ir in areas with lower resistivity values. In detail, the minimum resistivity and infiltration values correspond to point A5 and point A7 (50-70 Ωm and less than 0.05 m/day, respectively) and the maximum values with the point A3 (more than 500 Ωm and 3 m/day). The coefficient of determination (R 2 ) is below 0.9, but the positive slope of the relationship suggests that areas with ERT high values will have a high vertical hydraulic conductivity (k).

Assessment of Recharge Potential
The results of the 23 ERT cross-sections used for the characterization of the subsoil below the Ca n'Albareda meander ( Figure 1) show a high contrast in the electrical resistivity values, ranging from 5 to 1300 Ωm. This fact demonstrates the variability of the lithology and the hydraulic properties of sediments existing in the subsurface. Nevertheless, there is an overall trend of an increase in the electrical resistivity values from the bottom to the top and a noticeable variation in the upper levels.
By examining each ERT cross-section one by one, a layer with relatively low electrical resistivity values can be identified between the surface and about seven meters depth, and an intermediate layer with higher values. At the bottom, a low resistivity layer is identified, and it can be considered the impervious basement of the alluvial aquifer. The depth of this limit in the ERT cross-sections ranges between 15 (Rp1) and almost 40 m (Rp7) and, with the meander being a topographically flat area, defines an aquifer with a very variable thickness depending on the part of the meander considered.
In addition, as can also be seen at the bottom of Figure 9, the electrical resistivity response of the materials forming the aquifer is different due to the heterogeneity of alluvial sediments. Data from all ERT cross-sections have been interpolated to provide an indirect and continuous distribution of the hydraulic conductivity from the electrical resistivity values and the thickness all over the studied area. The kriging method has demonstrated to be well-adapted to the interpolation for electrical resistivity data [60], using the information provided by the variogram to better define aquifer/aquiclude boundaries and as an input in the kriging algorithm for interpolation throughout the studied area ( Figure 10). Data from all ERT cross-sections have been interpolated to provide an indirect and continuous distribution of the hydraulic conductivity from the electrical resistivity values and the thickness all over the studied area. The kriging method has demonstrated to be well-adapted to the interpolation for electrical resistivity data [60], using the information provided by the variogram to better define aquifer/aquiclude boundaries and as an input in the kriging algorithm for interpolation throughout the studied area ( Figure 10).
Once the geometry of the aquifer was obtained, the T parameter of Dar Zarrouk was calculated by multiplying the average electrical resistivity values with the thickness of the aquifer sediments from the water table to the upper limit of the aquiclude formation. The same interpolation algorithm (kriging) has also been used to obtain the distribution of the T parameter along the Ca n'Albareda meander. According to the geographical distribution of Dar Zarrouk's T parameter, a semiquantitative view of aquifer recharge potential has been obtained. The areas with the greatest potential for MAR are located in the east and south of the central area of the Ca n'Albareda meander and the areas with the least potential are located in the north of the meander (Figure 11). Once the geometry of the aquifer was obtained, the T parameter of Dar Zarrouk was calculated by multiplying the average electrical resistivity values with the thickness of the aquifer sediments from the water table to the upper limit of the aquiclude formation. The same interpolation algorithm (kriging) has also been used to obtain the distribution of the T parameter along the Ca n'Albareda meander. According to the geographical distribution of Dar Zarrouk's T parameter, a semiquantitative view of aquifer recharge potential has been obtained. The areas with the greatest potential for MAR are located in the east and south of the central area of the Ca n'Albareda meander and the areas with the least potential are located in the north of the meander ( Figure 11).   Once the geometry of the aquifer was obtained, the T parameter of Dar Zarrouk was calculated by multiplying the average electrical resistivity values with the thickness of the aquifer sediments from the water table to the upper limit of the aquiclude formation. The same interpolation algorithm (kriging) has also been used to obtain the distribution of the T parameter along the Ca n'Albareda meander. According to the geographical distribution of Dar Zarrouk's T parameter, a semiquantitative view of aquifer recharge potential has been obtained. The areas with the greatest potential for MAR are located in the east and south of the central area of the Ca n'Albareda meander and the areas with the least potential are located in the north of the meander (Figure 11).

Recharge Pond Management
After Ca n'Albareda's pond commissioning, problems due to low infiltration rate operation were identified and the facility came to a standstill. In surface MAR ponds-low-cost facilities by definition-the surface footprint and the cost of cleaning the basin must be minimized to maintain a good cost/profit ratio during its lifecycle [61]. We used the ERT technique to design the location and dimension of a pilot excavation, with the scope of minimizing basin rehabilitation costs through defining high-productivity areas [62]. The aim of the excavation at the pilot test area was to remove the shallowest layer formed by gravels with silt and clay matrix in the, a priori, most favorable area for artificial recharge. To eliminate the surface layer, it was necessary to carry out an excavation of approximately 0.8 m in depth. For the location of the pilot, the geophysical information of the NSZ characterization was critical since it was the only one available with a representative distribution of the heterogeneity of the subsoil of the infiltration pond. The sector was selected for having the highest electrical resistivity values in the NSZ campaign ( Figure 12) and the infiltration procedure could be restarted again.
After Ca n'Albareda's pond commissioning, problems due to low infiltration rate operation were identified and the facility came to a standstill. In surface MAR ponds-low-cost facilities by definition-the surface footprint and the cost of cleaning the basin must be minimized to maintain a good cost/profit ratio during its lifecycle [61]. We used the ERT technique to design the location and dimension of a pilot excavation, with the scope of minimizing basin rehabilitation costs through defining high-productivity areas [62]. The aim of the excavation at the pilot test area was to remove the shallowest layer formed by gravels with silt and clay matrix in the, a priori, most favorable area for artificial recharge. To eliminate the surface layer, it was necessary to carry out an excavation of approximately 0.8 m in depth. For the location of the pilot, the geophysical information of the NSZ characterization was critical since it was the only one available with a representative distribution of the heterogeneity of the subsoil of the infiltration pond. The sector was selected for having the highest electrical resistivity values in the NSZ campaign ( Figure 12) and the infiltration procedure could be restarted again.

Aquifer Recharge Potential Areas
The most favorable areas for artificial recharge must be able to store water and transmit it since this same water must be able to be captured downstream. The parameter of interest is hydraulic conductivity. The positive slope of the relationship between Ir and resistivity suggested that areas of continuous high electrical resistivity values might act as zones of preferential groundwater flow within the aquifer under suitable hydrologic conditionals. However, after ERT cross-sections were acquired, there are two key aspects to have been considered before moving on to the hydrogeological interpretation stage. The first is the well-known problem of equivalence-not uniqueness in the solution of the inverse problem-inherent in the interpretation of electrical resistivity over horizontally layered media, and the second is that ERT cross-sections are a 2D simplified representation of the subsurface.
Each of the ERT cross-sections used in the assessment of aquifer recharge has been taken to obtain an initial model and be able to address the equivalence problem. The midpoint apparent

Aquifer Recharge Potential Areas
The most favorable areas for artificial recharge must be able to store water and transmit it since this same water must be able to be captured downstream. The parameter of interest is hydraulic conductivity. The positive slope of the relationship between Ir and resistivity suggested that areas of continuous high electrical resistivity values might act as zones of preferential groundwater flow within the aquifer under suitable hydrologic conditionals. However, after ERT cross-sections were acquired, there are two key aspects to have been considered before moving on to the hydrogeological interpretation stage. The first is the well-known problem of equivalence-not uniqueness in the solution of the inverse problem-inherent in the interpretation of electrical resistivity over horizontally layered media, and the second is that ERT cross-sections are a 2D simplified representation of the subsurface.
Each of the ERT cross-sections used in the assessment of aquifer recharge has been taken to obtain an initial model and be able to address the equivalence problem. The midpoint apparent electrical resistivity values (one-dimensional (1D) data) have been inverted and the equivalence analysis has been performed using the RESIX-PLUS v2 software [63].
The results conform to three-layer models that describe K-type curves and are consistent with the information provided by the ERT cross-section in the Ca n'Albareda meander. K-type curves are characterized by a high electrical resistivity layer located between two more electrically conductive layers and can be unambiguously from the T Dar Zarrouk parameter.
The 1D high electrical resistivity layer is interpreted as the upper aquifer of Ca n'Albareda meander and the equivalence analysis has been performed by varying the h and ρ of the second layer within certain limits (Figure 13). Limits are marked by the error value of the inversion process, i.e., equivalent models with a similar error are acceptable. The relationship between h and ρ can be adjusted to a linear or potential function ( Figure 14) and the equivalent models of layer 2 have very similar Dar Zarrouk T values. the information provided by the ERT cross-section in the Ca n'Albareda meander. K-type curves are characterized by a high electrical resistivity layer located between two more electrically conductive layers and can be unambiguously from the T Dar Zarrouk parameter.
The 1D high electrical resistivity layer is interpreted as the upper aquifer of Ca n'Albareda meander and the equivalence analysis has been performed by varying the h and ρ of the second layer within certain limits ( Figure 13). Limits are marked by the error value of the inversion process, i.e., equivalent models with a similar error are acceptable. The relationship between h and ρ can be adjusted to a linear or potential function ( Figure 14) and the equivalent models of layer 2 have very similar Dar Zarrouk T values.   the information provided by the ERT cross-section in the Ca n'Albareda meander. K-type curves are characterized by a high electrical resistivity layer located between two more electrically conductive layers and can be unambiguously from the T Dar Zarrouk parameter. The 1D high electrical resistivity layer is interpreted as the upper aquifer of Ca n'Albareda meander and the equivalence analysis has been performed by varying the h and ρ of the second layer within certain limits ( Figure 13). Limits are marked by the error value of the inversion process, i.e., equivalent models with a similar error are acceptable. The relationship between h and ρ can be adjusted to a linear or potential function ( Figure 14) and the equivalent models of layer 2 have very similar Dar Zarrouk T values.   The standard deviation of the Dar Zarrouk T values is 4% lower than the average values. The results of the analysis validate the use of the Dar Zarrouk T parameter-already presented in the present paper-to characterize the aquifer unit of the Ca n'Albareda meander. Nevertheless, performing the procedure described has emphasized that a single determination of h and ρ would be difficult, if not impossible [64,65].

Conclusions
The ERT method has been a very useful tool for the characterization of the morphology of the alluvial aquifer of the Llobregat River, in particular for the characterization of surface recharge sites in the meander of Ca n'Albareda. This geophysical method has provided much more complete and representative information than direct techniques alone in an area of about 60 hectares. The geometry of the hydrogeological units and the aquifer-aquiclude contact, together with vertical and lateral variations, are accurately defined through the models derived from ERT cross-sections.
An empirical relationship has been derived from the electrical resistivity values and infiltration rates obtained at nine test sites, with an R 2 coefficient close to 0.9. The minimum resistivity and infiltration values are 50-70 Ωm and less than 0.05 m/day, respectively. On the other hand, the maximum values obtained are higher than 500 Ωm and 3 m/day. The positive slope of the relationship between infiltration and electrical resistivity suggested that areas of continuous high resistivity may act as zones of preferential flow within the aquifer under suitable hydrological conditions.
The intrinsic limitations of the methodology, such as the equivalence problem and the interpolation method used, can be solved by a simple analysis of the equivalent models and with direct information, obtained, for example, with lithological data from boreholes and hydraulic conductivity values estimated from pumping and/or infiltration tests.
Prior to the construction of recharge basins, it is highly recommended to conduct studies such as that presented in this paper to obtain a semiquantitative view at the catchment scale. From the point of view of water management, it is of particular importance to characterize variations in gravel layers. Free silty-clay matrix layers are related to areas of highest permeability, which are, therefore, the most suitable for aquifer artificial recharge procedures.