A Multidisciplinary Approach for Clarifying the Recharge Processes and Origin of Saline Water in the Semi-Arid Punata Alluvial fan in Bolivia

: The analysis of stable isotopes assisted in identifying that groundwater in the Punata alluvial fan is mainly recharged by heavy ﬂash ﬂoods, and the recharge from rainfall is of less importance. In addition, the hydrochemical analysis identiﬁed the Pucara River as the main source of recharge. Other streams in the north and northwest of the fan do not seem to contribute to the recharge. The hydrochemistry also shows that there is an increase of the Na + and Cl − concentrations in the middle and distal part of the fan. The salinization of groundwater is most likely a result of the mixing of fresh water with residual saline pore water in the lacustrine deposits and/or ion exchange within these layers. Geophysical surveys assisted in describing the aquifer system layering, and indicated a ﬁne-grained bottom layer where ion exchange might occur. This study demonstrates that the integration of several methods (e.g., hydrochemistry, hydrogeophysics, and stable isotopes) is valuable for clarifying ambiguities during the interpretation process and for characterizing hydrogeological processes in alluvial fans in general.


Introduction
Surface water access is considered to play a crucial role for the development of economic and social activities, but groundwater reserves can be as important as surface water and can thus play a similar role especially in arid and semiarid regions of the world [1][2][3]. Alluvial fans are generally located in arid and semi-arid zones, and they have the nature of having permeable layers where water can be stored, thus, they can be important sources of water supply. For sustainable management of groundwater, however, knowledge of origin, recharge process, flow direction, and mineralization is needed.
In Bolivia, there is generally a lack of comprehensive research which can provide knowledge of importance for the planning of further sustainable strategies for groundwater exploitation. The Punata where the question mark means that the depth to bedrock is still unknown, from Gonzales Amaya, et al. [15].
The regional geology is described by UNDP-GEOBOL [12], GEOBOL [16], and May, et al. [17]. The basement of the study area is formed of Ordovician rocks. These rocks are mainly limonite and where the question mark means that the depth to bedrock is still unknown, from Gonzales Amaya, et al. [15].
The regional geology is described by UNDP-GEOBOL [12], GEOBOL [16], and May, et al. [17]. The basement of the study area is formed of Ordovician rocks. These rocks are mainly limonite and quartz sandstone. Overlying Silurian formations consist of alternating strata of micaceous minerals and reddish and blackish shales. Dolostone, marlstone and limestone are the typical rocks in the Cretaceous formations, which overlie the Silurian deposits. The valleys are filled with unconsolidated sediments from the Quaternary, the main units are: fluvial, lacustrine, fluviolacustrine, and colluvial-alluvial sediments. These units form terraces on riverbanks and also beds and alluvial fans at the mouths of creeks and rivers. During the lower Pliocene, due to the uplift of the mountain massif, tectonic valleys and closed lakes were created. The climate in this period was predominately dry and saline material with high clay content was deposited [12]. Figure 1b shows the deposition of Quaternary sediments in the Punata alluvial fan, as follows: (1) In the apex of the fan, fluvial (Qa) and terrace sediments were deposited. These deposits are mainly composed of coarse material, varying from boulders to sand; (2) The middle of the fan is a transitory zone from alluvial (Qaa) to fluvio-lacustrine (Qfl) deposits. The grain size of these deposits ranges from gravel to silt and clay; (3) The distal part of the fan is mainly composed of fluvio-lacustrine deposits (Qfl) and the deposits are finegrained and mainly composed of silt and clay.

Water and Soil Sample Collection
This study included a total of 57 samples from two different sources: groundwater and surface water (Tables 1 and 2). The samples from groundwater are located within the Punata fan, while the latter are distributed along surface water bodies in the neighboring basins. The sampling depth of groundwater depended on the position of the well screens, which are often located at the bottom of the borehole. However, there are some wells that have screens at different depths. Table 2 provides information about the depth of sampling in wells and the aquifer unit where the bottom screen is located. The samples were spatially distributed in order to take into consideration all the possible recharge sources. At each sampling location, Electrical Conductivity (EC) and pH values were measured. The water samples were analyzed at the Center of Water and Environmental Sanitation (Centro de Aguas y Saneamiento Ambiental) within the Universidad Mayor de San Simon (UMSS) in Bolivia. The analyses include total dissolved solids (TDS), major ions such as calcium (Ca 2+ ), magnesium (Mg 2+ ), sodium (Na + ) bicarbonate (HCO 3 − ), chloride (Cl − ), sulfate (SO 4 2− ); and, minor species, such as potassium (K + ), and nitrate (NO 3 − ). The concentrations are given in milligram per liter (mg/L), while EC is presented in micro Siemens per centimeter (µS/cm). The water sampling followed the protocol and standards methods of APHA [18]. After purging the boreholes, the water samples were collected in 500 mL sterilized polyethylene bottles and stored below 4 • C until the analysis was performed. Bicarbonate was determined by H 2 SO 4 titration (with a mixed indicator bromocresol green-methyl red, and corroborated with an end-point pH measurement A total of 20 top soil samples were collected within the Punata fan ( Figure 2). The samples are mainly located in the head of the fan due to the fact that in this area it is expected to have a greater variation in the soil physical properties when compared to the middle and distal parts of the fan where the top soil is mostly clayey. In each location, soil samples were collected and infiltration tests were performed, both at three different depths: 30, 50 and 100 cm below the surface. The following physical properties were determined in the laboratory: texture, bulk density ( ), field capacity (θ fc ), wilting point (θ wp ), total porosity (n), and content of gravel, sand, clay, and silt. The infiltration surveys were performed with a Guelph permeameter, and when conditions were not suitable (i.e., coarse soil) the double ring infiltrometer was used. The estimation of the hydraulic conductivity (K) was done by applying the approaches that are defined in Reynolds and Elrick [20] and Elrick, et al. [21]. In order to verify the consistence of the estimates, these were compared with the values from the ROSETTA software [22]. ROSETTA uses a neural network analysis to estimate hydraulic parameter, such as K. This software can use limited information (i.e., grain size distribution) or more detailed information (i.e., grain size distribution, bulk density, porosity, wilting point, and/or field capacity) for the estimation of K.

Ion Chemistry and Chemical Weathering Process
Multiple techniques were used to visualize and analyze the hydrochemical processes occurring within the study area. Scatter plots and Chadha diagrams [23] were used to help in the task of distinguishing patterns and relationships in the aquifer system. Many studies, e.g., Han, et al. [24], Edmunds, et al. [25], and Herczeg, et al. [26], demonstrated that the concentration of major ions is a function of the direction and distance of groundwater flow. Weathering processes are important when interpreting the ion concentrations in the water sample. By performing chemical analysis of groundwater samples, it is possible to track the mineral weathering occurred along a groundwater path within an aquifer system. Appelo and Postma [27], Marion [28], Drever [29], and Datta and Tyagi [30] summarize some of the typical known weathering reactions, and how these reactions are associated to the weathering of some minerals and rocks (e.g., halite dissolution).
The ability of ion exchange is an important phenomenon and determines the mobility of some chemical elements. Ion exchange sites are often found in clay and organic soil materials, and some authors have proposed the following ordering of cation exchangeability [6,27]: Na + > K + > Mg 2+ > Ca 2+ , the divalent ions tend to replace the monovalent ions, but it is reversible and at high activities the monovalent ions can replace the divalent ions. The relation of [(Ca 2+ + Mg 2+ )-(HCO 3 − + SO 4 2− )] and (Na − -Cl − ) were used by several authors (e.g., Bouzourra, et al. [31], García, et al. [32], and Jalali [33]), for determining the gain or loss of (Ca 2+ + Mg 2+ ) due to silicate weathering and/or ion exchange, while (Na + -Cl − ) shows the gain or the loss of Na + provided by evaporates.

Stable Isotopes
The application of stable isotope analysis can provide a better understanding of some hydrogeological processes. The environmental isotopes, such as oxygen-18 ( 18 O) and deuterium ( 2 H), are non-anthropogenic and natural. They can assist in the understanding of hydrochemical problems, for instance, in the identification of recharge sources and flow patterns, and also in the development of a water budget [34]. The occurrence of 18 O and 2 H in groundwater might have a meteoric origin, therefore there is a strong bond between these isotope concentrations and the values of precipitation reflected in the Global Meteoric Water Line (GMWL); this relationship is given by a slope of 8 and a 2 H excess of +10‰ [35]. Deviations from the GMWL can be interpreted as being caused by local or regional climatic condition, or by precipitation that occurred during warmer or colder climate than at present. The main effects that cause deviations from the GMWL are temperature, amount effect, evaporation, and altitude [36].

Stable Isotopes Composition
The stable isotopic ratio of the water samples ranged from −14.5‰ to −9.6‰ for δ 18 O and from −106.5‰ to −77.9‰ for δ 2 H (refer to Table 1). Rainfall records from the researches [9,[37][38][39] are used for comparison, since it was not possible to collect rainfall samples within the study area. Stimson, Frape, Drimmie, and Rudolph [9] collected five rainfall water samples in the northwestern nearby basin of Cochabamba. The Cochabamba basin is located 40 km away, and has similar climate and physiographic characteristics. The rainfall samples from Cochabamba were collected during the months of January to June, and have been used for defining a local meteoric water line (LMWL): δ 2 H = 7.8δ 18 O + 11. In the study area the groundwater samples showed less variation in the isotopic composition than samples from surface water, and have similar values than the ones that are found in [37]. Both surface water and groundwater samples are shifted to the right of the GMWL (refer to Figure 3) highlighting the evaporation processes within the samples. A local evaporation water line (LWEL) was estimated: δ 2 H = 5.6δ 18 O − 24.

General Ions Trends
The physicochemical results of the 57 water samples are listed in Table 2. It can be noted that there are some variations among the collected samples, especially for EC, Na + , Cl − and TDS. This is most likely due to the fact that groundwater samples in the distal part of the fan (e.g., northwestern part) have longer flow paths, and have dissolved more minerals, which increase the TDS, and consequently, the EC.
The concentration of the major elements ranges widely, for instance, concentrations of Cl − and Na + can range from 0. 12 Table 3 summarizes the laboratory results that were obtained from the top soil samples. The grain size distribution of the soil samples decreases in size from the apex towards the distal part of the fan. In the same way, the values of K tend to be greater within the head of the fan and smaller in the distal and mid of the fan. Figure 2a shows the texture classification of the soil samples, where three different points clusters can be distinguished. The first of them have higher percentages of sand and all of these samples are located within the apex and the head of the fan. While the other two clusters have higher percentages of silt and clay, and they are located mainly within the middle and distal part of the fan. The difference in percentages might be explained due to the fact that energy of Pucara River is greater close to the apex, therefore carries coarser material. While towards the middle and distal part of the fan the energy of the river decreases, and therefore transports only finer material such as clay and silt.  Figure 2b displays the values and distribution of K at three different depths (i.e., 30, 50 and 100 cm). In general, the values of K are smaller at depths of 30 cm, and greater at depths of 50 and 100 cm. Also, it is possible to suggest that K values are greater in the apex and head of the fan, than the ones on the distal and middle parts of the fan. The difference between K values are mainly due to the fact that coarser material is more likely in the apex and head of the fan, which increases the values of K. the distal and mid of the fan. Figure 2a shows the texture classification of the soil samples, where three different points clusters can be distinguished. The first of them have higher percentages of sand and all of these samples are located within the apex and the head of the fan. While the other two clusters have higher percentages of silt and clay, and they are located mainly within the middle and distal part of the fan. The difference in percentages might be explained due to the fact that energy of Pucara River is greater close to the apex, therefore carries coarser material. While towards the middle and distal part of the fan the energy of the river decreases, and therefore transports only finer material such as clay and silt.  Figure 2b displays the values and distribution of K at three different depths (i.e., 30, 50 and 100 cm). In general, the values of K are smaller at depths of 30 cm, and greater at depths of 50 and 100 cm. Also, it is possible to suggest that K values are greater in the apex and head of the fan, than the ones on the distal and middle parts of the fan. The difference between K values are mainly due to the fact that coarser material is more likely in the apex and head of the fan, which increases the values of K.

Groundwater Recharge
The composition of stable isotopes in the Punata alluvial fan is similar to those generally found in semi-arid and arid regions, which are dominated by evaporation signatures, e.g., The rainfall that takes place in the study area come from the Amazon basin (i.e., origin in the Atlantic ocean), therefore the depleted values in rainfall samples (i.e., −12.1‰ to −8.4‰ for δ 18 O and −79.8‰ to −58.8‰ for δ 2 H), as found by [9], for the neighboring Cochabamba basin are expected since the study area is a region further inland from the Atlantic ocean. The LMWL of the study area ( Figure 3) is similar to the one found by [40] in the northern region of Chile (δ 2 H = 7.7δ 18 O + 9.7), and the one found in [39] as obtained from samples in a transect from the Bolivian Amazon to the Andes range (δ 2 H = 8.0δ 18 O + 15.2). The slope of these studies are lower or equal than eight from the GMWL, which probably is explained by an amount effect: during heavy rainfalls the isotopic composition tends to decrease with time due to a continuous isotopic exchange within the raindrops and the air mass below the storm [41]. Groundwater and surface water samples in the study areas have an evaporation trend, which is similar to the one found in Oruro [42], at 150 km west from the Punata region. The LEWL from Oruro has a slope of 5, and since both Punata and Oruro areas receive storms primarily from the Amazon basin, and Punata is somewhat more humid than Oruro, the LEWL in Punata (Figure 3) should have a slope slightly higher than 5 (i.e., 5.6). The LMWL slopes from Punata and Oruro fit with the description that in arid and semi−arid regions the evaporation trends have slopes as low as 4 to 6, highlighting the recycling of evaporated water into cloud masses [43].

Groundwater Recharge
The composition of stable isotopes in the Punata alluvial fan is similar to those generally found in semi-arid and arid regions, which are dominated by evaporation signatures, e.g., The rainfall that takes place in the study area come from the Amazon basin (i.e., origin in the Atlantic ocean), therefore the depleted values in rainfall samples (i.e., −12.1‰ to −8.4‰ for δ 18 O and −79.8‰ to −58.8‰ for δ 2 H), as found by [9], for the neighboring Cochabamba basin are expected since the study area is a region further inland from the Atlantic ocean. The LMWL of the study area ( Figure 3) is similar to the one found by [40] in the northern region of Chile (δ 2 H = 7.7δ 18 O + 9.7), and the one found in [39] as obtained from samples in a transect from the Bolivian Amazon to the Andes range (δ 2 H = 8.0δ 18 O + 15.2). The slope of these studies are lower or equal than eight from the GMWL, which probably is explained by an amount effect: during heavy rainfalls the isotopic composition tends to decrease with time due to a continuous isotopic exchange within the raindrops and the air mass below the storm [41]. Groundwater and surface water samples in the study areas have an evaporation trend, which is similar to the one found in Oruro [42], at 150 km west from the Punata region. The LEWL from Oruro has a slope of 5, and since both Punata and Oruro areas receive storms primarily from the Amazon basin, and Punata is somewhat more humid than Oruro, the LEWL in Punata (Figure 3) should have a slope slightly higher than 5 (i.e., 5.6). The LMWL slopes from Punata and Oruro fit with the description that in arid and semi−arid regions the evaporation trends have slopes as low as 4 to 6, highlighting the recycling of evaporated water into cloud masses [43]. The surface waters in Figure 3 are clustered in three groups, which probably are due the fact that samples were collected in different seasons. Samples collected during the dry period would tend to be more isotopically enriched due to evaporation effects, while samples that were collected during the rainy period might tend to be more depleted probably due to the amount effect explained above. Two surface water samples were clearly more depleted than the rest of the samples; the locations of these samples are springs that are located in the neighbor basin (Paracaya streams, refer to Figure  1a). These springs may be the discharge areas of water infiltrated within the secondary porosity of the surrounding Ordovician rocks, hence the isotopic composition cannot be explained by evaporation because they are located within the GMWL highlighting a meteoric origin. Therefore, one explanation might be that these waters have been infiltrated during different climate conditions where meteoric waters were more isotopically depleted. However, more investigation is needed in this area to fully understand the local hydrological processes.
Groundwater samples showed less variation regarding isotopic composition, with an approximately variation of 2‰ and 13‰ for δ 18 O and δ 2 H, respectively. Groundwater samples have similar isotopic composition with samples of surface water collected during the rainy season. The The surface waters in Figure 3 are clustered in three groups, which probably are due the fact that samples were collected in different seasons. Samples collected during the dry period would tend to be more isotopically enriched due to evaporation effects, while samples that were collected during the rainy period might tend to be more depleted probably due to the amount effect explained above. Two surface water samples were clearly more depleted than the rest of the samples; the locations of these samples are springs that are located in the neighbor basin (Paracaya streams, refer to Figure 1a). These springs may be the discharge areas of water infiltrated within the secondary porosity of the surrounding Ordovician rocks, hence the isotopic composition cannot be explained by evaporation because they are located within the GMWL highlighting a meteoric origin. Therefore, one explanation might be that these waters have been infiltrated during different climate conditions where meteoric waters were more isotopically depleted. However, more investigation is needed in this area to fully understand the local hydrological processes.
Groundwater samples showed less variation regarding isotopic composition, with an approximately variation of 2‰ and 13‰ for δ 18 O and δ 2 H, respectively. Groundwater samples have similar isotopic composition with samples of surface water collected during the rainy season. The rainy season in the study area is characterized by heavy storms of short duration, which usually originates flash flood events. Since the heavy storms may be more isotopically depleted, the originated flash floods would be depleted as well, which is in line with the similar isotopic composition with groundwater samples and surface water samples that were collected during rainy season. On the other hand, some groundwater samples located within the middle and distal part of the fan are more isotopically enriched, which might be due to vertical infiltration of evaporated water. However, this pattern is not reproduced entirely in all of the samples located in the distal and middle part of the fan, suggesting that the infiltration from the local streams during heavy floods is the predominant recharge process within the aquifer system. Hence, it might be suggested that light rainfalls, which are more enriched in δ 18 O and δ 2 H, are of minor importance in the recharge process. Consequently, heavy storms, which lead to heavy flash floods events, are of major importance for the recharge process in the Punata alluvial fan.

Origin of Saline Water
The water chemistry in the study area generally evolves from Ca 2+ /Na + -HCO 3 − to Na + -Cl − water type. Figure 4a shows that surface water have predominately higher concentrations of HCO 3 − , the arrows in this figure resemble the flow directions: water from rivers and streams infiltrates into the aquifer system and flows towards the distal part of the fan. Surface water from Pucara and Pocoata River have similar composition (i.e., Ca 2+ -HCO 3 − type), while Paracaya and San Benito surface waters have lower ion concentration and are more Na+ enriched (e.g., Na + -HCO 3 − water type).
In the distal part of the fan, the water type is dominated by higher concentrations of Na + and Cl − . rainy season in the study area is characterized by heavy storms of short duration, which usually originates flash flood events. Since the heavy storms may be more isotopically depleted, the originated flash floods would be depleted as well, which is in line with the similar isotopic composition with groundwater samples and surface water samples that were collected during rainy season. On the other hand, some groundwater samples located within the middle and distal part of the fan are more isotopically enriched, which might be due to vertical infiltration of evaporated water. However, this pattern is not reproduced entirely in all of the samples located in the distal and middle part of the fan, suggesting that the infiltration from the local streams during heavy floods is the predominant recharge process within the aquifer system. Hence, it might be suggested that light rainfalls, which are more enriched in δ 18 O and δ 2 H, are of minor importance in the recharge process. Consequently, heavy storms, which lead to heavy flash floods events, are of major importance for the recharge process in the Punata alluvial fan.

Origin of Saline Water
The water chemistry in the study area generally evolves from Ca 2+ /Na + -HCO3 − to Na + -Cl − water type. Figure 4a shows that surface water have predominately higher concentrations of HCO3 − , the arrows in this figure resemble the flow directions: water from rivers and streams infiltrates into the aquifer system and flows towards the distal part of the fan. Surface water from Pucara and Pocoata River have similar composition (i.e., Ca 2+ -HCO3 − type), while Paracaya and San Benito surface waters have lower ion concentration and are more Na+ enriched (e.g., Na + -HCO3 − water type). In the distal part of the fan, the water type is dominated by higher concentrations of Na + and Cl − . Four possible hypothesis might explain the increased salinization of the groundwater: (1) dissolution of halite; (2) infiltration of irrigation waters and/or excess on the use of fertilizers; (3) vertical diffuse recharge through the unsaturated soil; (4) or due to the existence of saline deposits lying in the lacustrine bottom of the local aquifer system. When hypothesis (1) was analyzed it would have been expected that the dissolution of halite in water released equal concentration of Na + and Cl − , but the results from this study deviates slightly from the expected 1:1 relation (Figure 4b). Most of the points are above the 1:1 line, which means that the concentration of Na + is higher than Cl − . Figure 4b also displays the saturation indexes (SI) with respect to halite, which are, on average, at least a factor of 8 less than the saturation point. Moreover, the local geology does not describe any halite deposit close to the study area. These indicators conduct to reject the hypothesis (1) that the dissolution of halite is the source of groundwater salinization.
For hypothesis (2) and (3); first, a spatial analysis of the ion concentrations in the study areas was done in order to identify the possible recharge areas and chemical evolution of groundwater samples. Figure 5 shows the circular diagrams with the concentration of the major ions in water When hypothesis (1) was analyzed it would have been expected that the dissolution of halite in water released equal concentration of Na + and Cl − , but the results from this study deviates slightly from the expected 1:1 relation (Figure 4b). Most of the points are above the 1:1 line, which means that the concentration of Na + is higher than Cl − . Figure 4b also displays the saturation indexes (SI) with respect to halite, which are, on average, at least a factor of 8 less than the saturation point. Moreover, the local geology does not describe any halite deposit close to the study area. These indicators conduct to reject the hypothesis (1) that the dissolution of halite is the source of groundwater salinization.
For hypothesis (2) and (3); first, a spatial analysis of the ion concentrations in the study areas was done in order to identify the possible recharge areas and chemical evolution of groundwater samples. Figure 5 shows the circular diagrams with the concentration of the major ions in water samples. The samples that were located in the northernmost part of the Punata fan are not displayed because surface water in overall has a similar chemical composition. Groundwater samples that were located in the discharge areas of Pucara and Pocoata Rivers have similar chemical composition as the surface samples; hence, this strengthens the isotopic analyses that suggest that the Punata groundwater is mainly recharged by river flows. Groundwater samples tend to increase the TDS towards the distal part of the fan. However, the northwestern border of the fan yielded the highest TDS, where Na + and Cl − are the dominant ions. If vertical infiltration would have been even in all of the study area, and since the crops and irrigation practices are, in general, the same in the Punata fan, it would be expected to find groundwater samples with more or less the same ion concentrations. Moreover, the most used method for irrigation in the agricultural lands of the study area is by flooding the terrains. Due to the semiarid conditions of the region, it would be expected that part of the irrigation water would be evaporated and/or transpired by plants, therefore salt ions would be deposited in the soil. After heavy rainfalls or when irrigation water is not evaporated nor transpired, it would be expected that water would infiltrate downwards, and hence, transporting all the solutes. However, the groundwater samples yielded very low concentrations with respect to NO 3 − (refer to Table 2), which indicates that there is no contamination from human activities, such as septic pits or fertilizers, even though such activities take place within the Punata fan. Therefore, the low concentration values of NO 3 − suggests that reduction processes occurred during infiltration, or that the water infiltrating vertically does not reach the sampled groundwater. Since the genesis of an alluvial fan entails an interfingering of different deposits with different grain sizes, it is probable that lenses with high concentration of clay in the unsaturated zone might act as a semipermeable barrier that prevents the vertical infiltration to reach the groundwater. Hence, it seems that hypothesis (2) and (3) do not explain the groundwater salinization in the study area. samples. The samples that were located in the northernmost part of the Punata fan are not displayed because surface water in overall has a similar chemical composition. Groundwater samples that were located in the discharge areas of Pucara and Pocoata Rivers have similar chemical composition as the surface samples; hence, this strengthens the isotopic analyses that suggest that the Punata groundwater is mainly recharged by river flows. Groundwater samples tend to increase the TDS towards the distal part of the fan. However, the northwestern border of the fan yielded the highest TDS, where Na + and Cl − are the dominant ions. If vertical infiltration would have been even in all of the study area, and since the crops and irrigation practices are, in general, the same in the Punata fan, it would be expected to find groundwater samples with more or less the same ion concentrations. Moreover, the most used method for irrigation in the agricultural lands of the study area is by flooding the terrains. Due to the semiarid conditions of the region, it would be expected that part of the irrigation water would be evaporated and/or transpired by plants, therefore salt ions would be deposited in the soil. After heavy rainfalls or when irrigation water is not evaporated nor transpired, it would be expected that water would infiltrate downwards, and hence, transporting all the solutes. However, the groundwater samples yielded very low concentrations with respect to NO3 − (refer to Table 2), which indicates that there is no contamination from human activities, such as septic pits or fertilizers, even though such activities take place within the Punata fan. Therefore, the low concentration values of NO3 − suggests that reduction processes occurred during infiltration, or that the water infiltrating vertically does not reach the sampled groundwater. Since the genesis of an alluvial fan entails an interfingering of different deposits with different grain sizes, it is probable that lenses with high concentration of clay in the unsaturated zone might act as a semipermeable barrier that prevents the vertical infiltration to reach the groundwater. Hence, it seems that hypothesis (2) and (3) do not explain the groundwater salinization in the study area. Another approach was used to assess hypothesis (2) and (3), which consisted in performing onedimensional (1D) models of water infiltration in the software Hydrus-1D [44]. Daily records of precipitation from 1 August 2010 to 31 July 2017 were used to simulate the water inflow in the 1D model. Hydrus-1D besides meteorological information (e.g., precipitation) required information Another approach was used to assess hypothesis (2) and (3), which consisted in performing one-dimensional (1D) models of water infiltration in the software Hydrus-1D [44]. Daily records of precipitation from 1 August 2010 to 31 July 2017 were used to simulate the water inflow in the 1D model. Hydrus-1D besides meteorological information (e.g., precipitation) required information about soil properties. The soil samples showed that it is very likely to find a semipermeable top soil up to two meters thick composed mainly of clay and sand, with K ranging from 15 to 206 cm/day. On the other hand, pumping tests yielded K between 100 and 900 cm/day, which corresponds to layers mainly composed of coarse material [12]. Figure 6a depicts the simulation of the water content in a soil column, where it was assumed that there is a semipermeable top layer with a K of 60 cm/days (average of K values from the soil samples), the thickness of this top layer is 2 m. The rest of the column is assumed to have a K value of 500 cm/day (average K from the research of UNDP-GEOBOL [12]). Figure 6b has the same soil column characteristic as in Figure 6a, with the addition of a semipermeable layer at a depth of 10 m. This additional layer has the same properties as the top layer. about soil properties. The soil samples showed that it is very likely to find a semipermeable top soil up to two meters thick composed mainly of clay and sand, with K ranging from 15 to 206 cm/day. On the other hand, pumping tests yielded K between 100 and 900 cm/day, which corresponds to layers mainly composed of coarse material [12]. Figure 6a depicts the simulation of the water content in a soil column, where it was assumed that there is a semipermeable top layer with a K of 60 cm/days (average of K values from the soil samples), the thickness of this top layer is 2 m. The rest of the column is assumed to have a K value of 500 cm/day (average K from the research of UNDP-GEOBOL [12]). Figure 6b has the same soil column characteristic as in Figure 6a, with the addition of a semipermeable layer at a depth of 10 m. This additional layer has the same properties as the top layer.  Figure 1a for location). Type of climate and rate of sedimentation are also displayed. After UNDP−GEOBOL [12]. Figure 6a,b shows that at the surface (i.e., 0 m) there are great variation in the soil moisture: after each event of precipitation, there is an increase in the soil water content, however during the dry season and due to evaporation effects the soil water content reaches the water residual content point. At a depth of 30 m, the soil water content, after four years of precipitation, barely reaches the value of 0.1, which means that the soil water content has not reached the saturation point. At the bottom of the profiles, the soil moisture remains constantly within the residual water content. The small values of the soil water content between 30 and 50 m depth highlights the assumption that water from rainfall/irrigation is not a major process within the groundwater salinization and recharge process. Therefore, it is suggested that the hypotheses (2) infiltration of irrigation waters and/or excess use of fertilizers, and (3) vertical diffuse recharge through the unsaturated soil are of minor importance in the groundwater salinization and vertical recharge. This is in line with the isotopic analysis where it is suggested that rainfall has a low role in the groundwater recharge and infiltration from rivers are the main recharge sources.
The final proposed hypothesis that could explain the high concentration of Na + and Cl − ions is hypothesis (4) that suggests the existence of saline deposits, which led to a mixing of groundwater with residual pore water in the lacustrine clays and ion exchange. Figure 6c displays the lithology of borehole P305, which is located within the head of the fan (refer to Figure 1a). The lithology description shows a decrease of grain size with depth, and a thick layer at the bottom of high clay content. In the study area, during the early Pleistocene the climate was mainly dry, therefore successive events of evaporation occurred within the former paleolakes leading to a slow sedimentation where fine particles were deposited. Also, the conditions during that period might  Figure 1a for location). Type of climate and rate of sedimentation are also displayed. After UNDP−GEOBOL [12]. Figure 6a,b shows that at the surface (i.e., 0 m) there are great variation in the soil moisture: after each event of precipitation, there is an increase in the soil water content, however during the dry season and due to evaporation effects the soil water content reaches the water residual content point. At a depth of 30 m, the soil water content, after four years of precipitation, barely reaches the value of 0.1, which means that the soil water content has not reached the saturation point. At the bottom of the profiles, the soil moisture remains constantly within the residual water content. The small values of the soil water content between 30 and 50 m depth highlights the assumption that water from rainfall/irrigation is not a major process within the groundwater salinization and recharge process. Therefore, it is suggested that the hypotheses (2) infiltration of irrigation waters and/or excess use of fertilizers, and (3) vertical diffuse recharge through the unsaturated soil are of minor importance in the groundwater salinization and vertical recharge. This is in line with the isotopic analysis where it is suggested that rainfall has a low role in the groundwater recharge and infiltration from rivers are the main recharge sources.
The final proposed hypothesis that could explain the high concentration of Na + and Cl − ions is hypothesis (4) that suggests the existence of saline deposits, which led to a mixing of groundwater with residual pore water in the lacustrine clays and ion exchange. Figure 6c displays the lithology of borehole P305, which is located within the head of the fan (refer to Figure 1a). The lithology description shows a decrease of grain size with depth, and a thick layer at the bottom of high clay content. In the study area, during the early Pleistocene the climate was mainly dry, therefore successive events of evaporation occurred within the former paleolakes leading to a slow sedimentation where fine particles were deposited. Also, the conditions during that period might have led to evaporated deposits or residual pore water in the bottom of the aquifer system, but when freshwater entered in contact with the evaporated deposits the first ions to be mobilized were Na + and Cl − . Therefore, both ion exchange and the mixing of residual pore water in the lacustrine layers with fresh groundwater might explain the high concentration of Na + and Cl − .
In order to confirm and assist in the interpretation of hypothesis (4), several geophysical campaigns were performed. The methods used were Electrical Resistivity Tomography (ERT) and Transient Electromagnetics (TEM), for more details about the surveys layout and results refer to Gonzales Amaya, Dahlin, Barmen, and Rosberg [15] and Gonzales Amaya, et al. [45]. Figure 7 shows the electrical resistivity distribution and the corresponding hydrogeological interpretation from two different areas in the Punata fan: area 1 (A-1) and area 2 (A-2), refer to Figure 1a for location. A-1 is located in the eastern part of the fan within the head and middle part of the fan, while A-2 is located in the northwestern part of the fan within the middle and distal part of the fan. The interpretation of the geophysical results showed a thick layer of coarse material within the head and middle part of the fan, while towards the distal part of the fan, the materials are predominately finer. Within the higher resistivity top layers, it is possible to observe areas of lower resistivity, which are probably thin lenses of sand and clay. As explained before, these layers might act as semipermeable or impermeable barriers preventing the vertical water infiltration (i.e., vertical transport of solutes). Moreover, both areas highlighted a very low resistivity layer close to the bottom of the aquifer system. The low resistivity values (as low as 1 Ωm in A-1, and 0.1 Ωm in A-2) might be explained by saline water stored within the finer sediments and also by the residual pore water within the lacustrine clay. Hence, the location of this saline layer indicates that hypothesis (4) can explain the origin of saline groundwater reported in some wells within the Punata alluvial fan. have led to evaporated deposits or residual pore water in the bottom of the aquifer system, but when freshwater entered in contact with the evaporated deposits the first ions to be mobilized were Na + and Cl − . Therefore, both ion exchange and the mixing of residual pore water in the lacustrine layers with fresh groundwater might explain the high concentration of Na + and Cl − . In order to confirm and assist in the interpretation of hypothesis (4), several geophysical campaigns were performed. The methods used were Electrical Resistivity Tomography (ERT) and Transient Electromagnetics (TEM), for more details about the surveys layout and results refer to Gonzales Amaya, Dahlin, Barmen, and Rosberg [15] and Gonzales Amaya, et al. [45]. Figure 7 shows the electrical resistivity distribution and the corresponding hydrogeological interpretation from two different areas in the Punata fan: area 1 (A-1) and area 2 (A-2), refer to Figure 1a for location. A-1 is located in the eastern part of the fan within the head and middle part of the fan, while A-2 is located in the northwestern part of the fan within the middle and distal part of the fan. The interpretation of the geophysical results showed a thick layer of coarse material within the head and middle part of the fan, while towards the distal part of the fan, the materials are predominately finer. Within the higher resistivity top layers, it is possible to observe areas of lower resistivity, which are probably thin lenses of sand and clay. As explained before, these layers might act as semipermeable or impermeable barriers preventing the vertical water infiltration (i.e., vertical transport of solutes). Moreover, both areas highlighted a very low resistivity layer close to the bottom of the aquifer system. The low resistivity values (as low as 1 Ωm in A-1, and 0.1 Ωm in A-2) might be explained by saline water stored within the finer sediments and also by the residual pore water within the lacustrine clay. Hence, the location of this saline layer indicates that hypothesis (4) can explain the origin of saline groundwater reported in some wells within the Punata alluvial fan. For summarizing all of the findings of this research, a hydrogeological conceptual model is displayed in Figure 8. The hydrochemical and isotopic composition of surface waters is similar to those groundwater samples that were found in the discharge areas of Pucara and Pocoata River, highlighting the hypothesis that river flow is the main recharge source. Moreover, samples of For summarizing all of the findings of this research, a hydrogeological conceptual model is displayed in Figure 8. The hydrochemical and isotopic composition of surface waters is similar to those groundwater samples that were found in the discharge areas of Pucara and Pocoata River, highlighting the hypothesis that river flow is the main recharge source. Moreover, samples of Paracaya and San Benito streams have different isotopic signature (i.e., more depleted) and ion composition (i.e., lower TDS) than groundwater samples, suggesting that these streams do not contribute to the Punata groundwater recharge. The high concentrations of Na + and Cl − of some samples indicate a salinization process within the aquifer system. At the aquifer system, bottom geophysical surveys highlighted a saline zone and a layer with high clay content. Therefore, this saline zone and high clay content layer explains the occurrence of ion exchange and that residual saline pore water are the main origins of groundwater salinization. Due to the semiarid conditions of the Punata region, evapotranspiration might had led to the deposition of solutes that could have been transported during vertical water infiltration. However, the possible location of impermeable or semipermeable lenses formed of clay and sandy clay hindered the infiltration of solutes into the groundwater. Finally, the different depths of boreholes in the northwestern part of the fan (e.g., P30 and P19 in Figure 8) can assist in explaining the higher values of TDS in some wells. For instance, P19 is deeper than P30, but the latter has higher TDS, which might be due to the fact that the P30 is closer to the aquifer bottom where the salinization activity is higher. On the other hand, P5, which is located closer to the recharge area of Pocoata River, has a lower value of TDS, and might be due to the fact that there is a thicker cap layer above the saturated zone and the bottom of the borehole is not close to the saline layer. Paracaya and San Benito streams have different isotopic signature (i.e., more depleted) and ion composition (i.e., lower TDS) than groundwater samples, suggesting that these streams do not contribute to the Punata groundwater recharge. The high concentrations of Na + and Cl − of some samples indicate a salinization process within the aquifer system. At the aquifer system, bottom geophysical surveys highlighted a saline zone and a layer with high clay content. Therefore, this saline zone and high clay content layer explains the occurrence of ion exchange and that residual saline pore water are the main origins of groundwater salinization. Due to the semiarid conditions of the Punata region, evapotranspiration might had led to the deposition of solutes that could have been transported during vertical water infiltration. However, the possible location of impermeable or semipermeable lenses formed of clay and sandy clay hindered the infiltration of solutes into the groundwater. Finally, the different depths of boreholes in the northwestern part of the fan (e.g., P30 and P19 in Figure 8) can assist in explaining the higher values of TDS in some wells. For instance, P19 is deeper than P30, but the latter has higher TDS, which might be due to the fact that the P30 is closer to the aquifer bottom where the salinization activity is higher. On the other hand, P5, which is located closer to the recharge area of Pocoata River, has a lower value of TDS, and might be due to the fact that there is a thicker cap layer above the saturated zone and the bottom of the borehole is not close to the saline layer.

Conclusions
This study used the integration of the hydrochemistry of major elements, analysis of stable isotopes, modeling of 1D soil water movement, and geophysical surveys for providing a more comprehensive understanding of the recharge process and saline evolution in the Punata alluvial fan aquifer. Concentrations of Na + and Cl − suggest an increase of salinization towards the distal part of the fan. Both saline residual pore water in the lacustrine sediments and ion exchange might cause groundwater to become richer with respect to Na + and Cl − , rather than halite dissolution and transport of solutes from the unsaturated zone. The composition of δ 18 O and δ 2 H in groundwater samples are less enriched than those from precipitation. This suggests that groundwater in the Punata fan is recharged mainly by heavy flash floods. The vertical infiltration is of less importance during the recharge process, which is probably due to semipermeable lenses within the unsaturated zone. The hydrochemical analysis also suggests that the Pucara and Pocoata rivers recharge the Punata aquifer, while the northern streams do not recharge the aquifer.

Conclusions
This study used the integration of the hydrochemistry of major elements, analysis of stable isotopes, modeling of 1D soil water movement, and geophysical surveys for providing a more comprehensive understanding of the recharge process and saline evolution in the Punata alluvial fan aquifer. Concentrations of Na + and Cl − suggest an increase of salinization towards the distal part of the fan. Both saline residual pore water in the lacustrine sediments and ion exchange might cause groundwater to become richer with respect to Na + and Cl − , rather than halite dissolution and transport of solutes from the unsaturated zone. The composition of δ 18 O and δ 2 H in groundwater samples are less enriched than those from precipitation. This suggests that groundwater in the Punata Water 2018, 10, 946 15 of 17 fan is recharged mainly by heavy flash floods. The vertical infiltration is of less importance during the recharge process, which is probably due to semipermeable lenses within the unsaturated zone. The hydrochemical analysis also suggests that the Pucara and Pocoata rivers recharge the Punata aquifer, while the northern streams do not recharge the aquifer.
The multidisciplinary approach used in this research showed good results, because when ambiguities occurred during the interpretation process the use of other approaches helped in the final interpretation. Therefore, this study highlighted the importance of using several methods for improving the interpretation when the existing data is limited. The results that were obtained in this study will contribute in the making of policies for sustainable groundwater management in the Punata alluvial fan. However, in order to improve the knowledge of the recharge process and to verify the presented conclusions, it is recommended to take samples of the heavy flash floods and rainfall that produced these floods, and conduct an analysis of δ 18 O and δ 2 H composition, which includes monthly weighted rainfall and river flow. Also, it is important to conduct a comprehensive study of the mineralogy of rocks and soils within the aquifer system and nearby basin to confirm the weathering processes and mineralization processes.