Hydrogeochemical characterization of an evaporite karst area affected by sinkholes ( Ebro Valley , NE Spain )

The main processes controlling the hydrochemistry of an alluvium-covered evaporite karst area with high sinkhole risk (Ebro Valley, NE Spain) are examined by means of multivariate analyses (Principal Component Analysis and Hierarchical Cluster Analysis), ion correlations and geochemical speciation-solubility calculations. The hydrogeochemistry of the studied system seems to be governed by the interaction between the groundwater from the salt-bearing evaporitic karst aquifer and from the overlying Ebro River alluvial aquifer. The observed hydrochemical features in the alluvial-karst aquifer system are mainly determined by the relative contribution of gypsum/anhydrite and halite dissolution, showing a wide spectrum from relatively fresh recharge waters (mainly irrigation waters) to highly evolved groundwater from the evaporitic aquifer. The variability of these contributions is especially evident at sinkhole ponds which, in some cases, seem to be associated with discharge areas of the karst aquifer in the valley bottom alluvium. Calculated saturation indexes suggest that, in contrast to gypsum, the amounts of halite in the sampled portions of evaporitic aquifer are not large enough to attain equilibrium, which is consistent with the predominance of gypsum/anhydrite reported for these materials. Furthermore, the observed Na:Cl and Ca:SO4 correlations and stoichiometries suggest that other possible processes, such as glauberite dissolution or Na/Ca-exchange, generally play a minor role (compared to halite and gypsum dissolution) in this system. Another important process in the system is the dissolution of carbonate minerals (dolomite and, possibly, calcite) fostered by the input of CO2(g), which is probably produced by pedogenic processes. Dolomite dissolution seems to be particularly relevant in the evaporitic materials probably due to dedolomitisation triggered by gypsum/anhydrite dissolution. Salt karst. Alluvial aquifer. Geochemical modeling. Hydrochemistry.


INTRODUCTION
The subsurface dissolution of soluble minerals in carbonate and evaporitic materials by groundwater may cause the gravitational deformation and internal erosion of the overlying sediments, eventually leading to the settlement of the ground surface.Commonly, the geomorphic expression of these hazardous subsidence phenomena corresponds to closed depressions designated as sinkholes or dolines (Williams, 2004;Beck, 2005;Waltham, 2005;Gutiérrez et al., 2008b).
The Ebro valley in the outskirts of Zaragoza city (NE Spain) is most probably the area in Europe where subsidence risk related to evaporite dissolution has the greatest economic impact.In this sector, sinkholes related to the karstification of the highly soluble Tertiary Zaragoza Formation (gypsum/anhydrite, halite, glauberite) underlying Quaternary alluvium show a high spatial and temporal frequency (Gutiérrez et al., 2008a;Galve et al., 2009).The occurrence and reactivation of sinkholes in this mantled karst setting have lead to the demolition of many buildings and affect numerous linear infrastructures, including the high-speed Madrid-Barcelona railway (Guerrero et al., 2008b;Galve et al., 2009;Gutiérrez et al., 2009).
The effective management of the subsidence risk requires a detailed knowledge of the processes and factors involved in sinkhole development.Among other features, the mineralogy, hydrochemistry and hydrogeology of the target areas need to be well-known in order to understand the processes leading to karstification (Lamont-Black et al., 2002).Nevertheless, the vast majority of earlier works dealing with these processes have been carried out in carbonate karst aquifers (e.g.López-Chicano et al., 2001;Aquilina et al., 2003Aquilina et al., , 2005Aquilina et al., , 2006 and references therein; Barbieri et al., 2005;Tuccimei et al., 2005;Moore et al., 2009;Barberá and Andreo, 2012 and references therein;Bicalho et al., 2012).Most of the studies on evaporite karst terrains affected by sinkholes are focused on geomorphological aspects and the associated environmental problems (e.g.Klimchouk et al., 1996;Calaforra and Pulido-Bosch, 1999;Johnson and Neal, 2003;Gutiérrez et al., 2008c;Cooper and Gutiérrez, 2013 and references in all of them), and works addressing the hydrochemical and geochemical factors controlling the dissolution processes involved in sinkhole development are quite scarce (e.g.Kaçaroglu et al., 2001;Günay, 2002;Land, 2003;Lamont-Black et al., 2005;Yechieli et al., 2006;Omelon et al., 2006;Chiesi et al., 2010;Fidelibus et al., 2011;Apaydin and Aktas, 2012).
In spite of the relevant engineering and economic implications of the sinkhole hazard related to evaporite dissolution in the Ebro Valley, the relationships between the hydrogeology and geochemistry of the area and the dissolution phenomena remain poorly understood (Sánchez et al., 2004;Jiménez-Torrecilla et al., 2004).Moreover, to the best of our knowledge, the hydrogeochemical features and controls of this evaporite karst have not been comprehensively addressed by any earlier investigation.Thus, the main goal of this study is to obtain a better understanding of the hydrogeology and hydrogeochemistry of this area in relation to the features and processes involved in the development of sinkholes.The hydrogeochemistry at several types of water points (springs, wells, sinkhole ponds and ditches) is thoroughly examined with the assistance of multivariate analyses, examination of ion correlations and geochemical speciation-solubility calculations.Moreover, the relationships between the observed hydrochemical features, possible water-rock interaction processes, the composition of the evaporitic bedrock, and hydrogeological variables, are analysed.As a result, a general hydrochemical model is proposed, which will constitute the basis for future detailed studies aimed at exploring how this type of investigations may contribute to a more efficient analysis and management of the sinkhole risk.

Geological setting and climate
The studied 60km long reach of the Ebro Valley is located in the central sector of the Ebro Tertiary Basin, which constitutes the southern foreland basin of the Pyrenees (Fig. 1A).Here, the Ebro Valley has been carved in subhorizontally lying evaporites of the Oligo-Miocene Zaragoza Formation, deposited in an extensive high-salinity playa-lake (Quirantes, 1978, Ortí andSalvany, 1997).This evaporitic formation reaches more than 850m in thickness (Torrescusa and Klimowitz, 1990).In the subsurface the formation is primarily composed of anhydrite (CaSO 4 )/gypsum (CaSO 4 •2H 2 O), halite (NaCl) and glauberite (Na 2 Ca [SO 4 ] 2 ), with interlayered marls and mudrocks including calcite (CaCO 3 ), dolomite (CaMg(CO 3 ) 2 ) and quartz (SiO 2 ) grains (Salvany et al., 2007 and references therein;Guerrero et al., 2008a,b, in press).Halite and glauberite beds situated close to the surface, reach as much as 75m and 30m in thickness, respectively (Salvany, 2009).In outcrops, the evaporitic succession exhibits around 300m of laminated and nodular secondary gypsum derived from the hydration of anhydrite or the replacement of glauberite (Fig. 1B).The evaporitic formation changes laterally into impervious and insoluble clay facies in the eastern sector of the studied area (Fig. 1A).
In the studied reach of the valley, the terrace and pediment deposits overlying the evaporitic materials may reach 80m and fill kilometer-scale basins generated by dissolution-induced syn-sedimentary subsidence (Guerrero et al., 2013).These alluvial deposits mainly consist of gravels dominated by carbonate and siliciclastic materials with a sand-silt matrix, frequently cemented by carbonates.
Geomorphological studies and the analysis of paleosubsidence structures reveal that sinkhole formation in this area has been related to: i) progressive sagging resulting from interstratal karstification of glauberite and halite units; ii) upward propagation of cavities developed within the evaporite bedrock by progressive rock collapse; iii) upward propagation of voids through the alluvial mantle by soil cavity collapse and downward migration of detrital sediments into dissolutional voids (Gutiérrez et al., 2008a, Guerrero et al., 2008a, b, 2013).Some of the sinkholes in the studied area are permanently filled with phreatic waters (Fig. 2) and, therefore, represent a valuable source of information on the hydrochemical and hydrogeological features of the alluvial-karst aquifer system.
The climate in this region is Mediterranean with strong continental influence, characterized by hot summers A) Location and geological settings of the study area and position of the selected sampling points for hydrochemical characterisation.B) Schematic cross-section perpendicular to the Ebro Valley downstream of Zaragoza city, showing the general distribution of the main lithostratigraphic units and the weathering zone in which anhydrite and glauberite are replaced by secondary gypsum.The section is based on borehole data from the southern margin of the valley (Salvany, 2009) and assummes a roughly symmetric configuration at both sides of the valley, excavated along the axial zone of a very gentle syncline.A c t a , 1 1 ( 4 ) , 3 8 9 -4 0 7 ( 2 0 1 3 )  D O I : 1 0 . 1 3 4 4 / 1 0 5 . 0 0 0 0 0 2 0 5 2 Hydrogeochemistry of evaporite Karst in the Ebro valley (NE Spain) and cold winters.The average annual precipitation in Zaragoza area is around 340mm and the mean monthly temperature ranges from 9ºC to 21ºC (Ninyerola et al., 2005).Precipitation shows a considerable inter-annual variability and a significant proportion corresponds to intense, short-duration convective storms that commonly occur in spring and autumn.The year 2011, corresponding to the sampling campaing carried out specifically for this study, showed a total precipitation of 315mm, indicating a slightly dryer year in relation to the average.The mean annual reference evapotranspiration in Ebro Valley exceeds 1150mm, indicating a negative water balance of more than 800mm and, therefore, a climate characterized by semiarid conditions (Diputación General de Aragón, 2007).

Hydrogeology of the studied area
The hydrogeology of the area is characterized by two main interconnected and unconfined aquifers: i) the alluvial aquifer, composed by the alluvium underlying the valley bottom and the associated terrace deposits at the valley margins; ii) the karst aquifer developed in the salt-bearing evaporitic bedrock, which discharges in the previous one.The Ebro River and the associated valley bottom alluvium act as the base level for the alluvial-karst aquifer system (Durán et al., 2005) (Fig. 1A, B).
Very little is known about the karst aquifer developed in the fractured and karstified evaporitic bedrock.Although some authors assumed that it behaves as an aquitard (e.g.Martín, 1993), multiple lines of evidence indicate that is has a considerable permeability related to dissolution processes (Gutiérrez et al., 2007;Guerrero et al., 2008b;Galve et al., 2009).The recharge of the evaporitic aquifer occurs through infiltration of irrigation and precipitation water in bedrock outcrops and in local perched alluvial deposits located on the valley flanks.The groundwater that flows through the saltbearing evaporitic bedrock, attaining a progressively higher mineralization, ultimately discharges into the floodplain alluvium, as evidenced by the frequent presence in this aquifer of highly-saline waters associated with springs, wells and sinkhole ponds.Some authors have proposed that there may be areas of preferential discharge of the evaporite karst aquifer into the overlying alluvial aquifer (Sánchez et al., 2004;Jiménez-Torrecilla et al., 2004;Gutiérrez et al., 2007).
The water balance for the Ebro alluvial aquifer in the studied reach of the valley indicates that around 90% of the  recharge is related to irrigation (Durán et al., 2005).This is the main factor controlling the seasonal variations in the watertable, which may rise several meters on average by the end of the irrigation season (August, September), during the low flow period of the Ebro river.

Selection of water points for hydrochemical characterisation
In order to characterise the main hydrochemical features of the Ebro alluvial aquifer and to explore its relationship with the evaporitic karst aquifer in the study area, time series of hydrochemical data from 25 sampling points were included in this study (Table 1).The selection of water points was intended to cover the whole study area and to represent different types of water points (wells, springs, sinkhole ponds, ditches and the Ebro River).As displayed in Table 1, the hydrochemical information from some of the points could be obtained from the public database of the Ebro Waters Authority (Confederación Hidrográfica del Ebro; CHE).For the rest of sampling points included in this work, between 3 and 7 sampling surveys were carried out between January 2011 and January of 2012.During these surveys, some additional samples were collected for analysis at several of the points from the CHE database in order to confirm and complete their reported hydrochemistry.
Main features of the selected water points sampled for the hydrochemical analyses of the study area.The location of the points is shown in Figure 1.The codes of the sampling points consist of two capital letters followed by a number.The first capital letter indicates if the point corresponds to a well (W--), a spring (S--) or a pond (H--).The second letter is used to indicate if the water at that point can be considered as relatively concentrated (-A-) or diluted (-B-) compared to the average TDS of the study area, which is around 1600 mg/l.Additionally, IR points correspond to irrigation waters.Finally, the numbers allow differentiating the different points of each type  A c t a , 1 1 ( 4 ) , 3 8 9 -4 0 7 ( 2 0 1 3 )  D O I : 1 0 . 1 3 4 4 / 1 0 5 . 0 0 0 0 0 2 0 5 2 Hydrogeochemistry of evaporite Karst in the Ebro valley (NE Spain)

G e o l o g i c a
394 As displayed in Table 1, two of the wells included in this study withdraw water from the evaporitic aquifer (WA2 and WA7), whereas the rest of them are shallower and withdraw groundwater from the alluvial aquifer, although some influence of the evaporitic groundwaters is also probable.In the case of springs and sinkhole ponds, their predominant association with one aquifer or the other is uncertain and will be assessed on the basis of their hydrochemistry.

Methods for sampling and analysis
In the sampled wells, water was continuously pumped while monitoring the variations of pH, conductivity and temperature until their stabilization (within measurement error) before measuring the field parameters and taking the water samples for analyses.In the case of the sinkhole ponds, sampling surveys carried out in April 2009 demonstrated no variations in the hydrochemistry associated with stratification in the water column or the distance to the center of the pond.Accordingly, only one sample per pond was taken to be representative enough of the hydrochemistry at each of these points.
At each sampling point, separated samples for anion and cation analyses were taken and stored at 4ºC until analyses.The elements included in the chemical analyses (Na, Ca, Mg, K, Cl, sulfate and alkalinity) are the major ions in the studied system, according to the available geological and mineralogical information.Based on duplicate analytical determinations, the uncertainty associated with the reported values for these elements in the samples specifically taken for this study is estimated to be generally within ±8%.Samples for cation analyses were filtered through 0.1mm and acidified to a pH lower than 1 with ultrapure HNO 3 , whereas samples for anion and alkalinity analyses were stored unfiltered and unacidified.Temperature was measured in situ with a thermometer (accuracy 0.1ºC).The pH value was determined by a combined glass electrode with automated temperature correction, calibrated on a regular basis with standard buffer solutions of pH 4, 7 and 10.The reported accuracy for pH measurements is ±0.05 pH units.Electrical conductivity was measured by using a conductivity meter with automated temperature correction.
Alkalinity and chloride contents were determined by titration with a Mettler titrator with end-point electrode.Alkalinity measurements were carried out less than 24 hours after sampling.According to the hydrogeochemical characteristics of our samples, total and carbonate alkalinity can be considered as equivalent.Sulfate concentrations were determined by colorimetry using a modification of the Nemeth method (Nemeth, 1963).Sulfate and chloride concentrations were analysed within 3-4 weeks after sample collection.Inductively Coupled Plasma-Atomic Emission Spectrometry (ICP-AES) was used for the analyses of cations (Ca, Mg, Na and K).
Geochemical calculations that will be discussed throughout this work have been performed with the PHREEQC code (Parkhust and Appelo, 2013) using the WATEQ4F thermodynamic database (Ball and Nordstrom, 2001) included in the PHREEQC package.Thermodynamic data for glauberite from the pitzer.datdatabase distributed with the PHREEQC code have been added to the WATEQ4F database.The calculated charge balance error for the reported samples, as calculated with the PHREEQC code, is generally below 5% and never above 10%.

Methods for multivariate statistical analysis
Two different multivariate statistical analysis techniques were applied to the study of the hydrochemical dataset generated during this work: Principal Component Analysis (PCA) and Hierarchical Cluster Analysis (HCA).These techniques have been extensively applied to different types of geological problems and, particularly, they have been used to identify processes or environmental features controlling the behaviour of large hydrochemical datasets in different types of geological systems, including karst systems (López-Chicano et al., 2001;Moore et al., 2009 and references therein;Dassi, 2011;Barbera and Andreo, 2012;Bicalho et al., 2012 and references therein; among many others).
For principal component analysis, all hydrochemical data (except pH) were log-transformed prior to the statistical analyses, so that they more closely correspond to normallydistributed populations, and they were also autoscaled by calculating their standard scores (z scores).For hierarchical cluster analysis, different similarity/dissimilarity measurements and linkage methods have been tested, and the results suggest that Euclidean distance as similarity measurement and Paired group linkage and Ward's method produced the most distinctive groups.Dendrograms have been interpreted based on visual examination and according to the water types described previously.
Both types of analysis have been applied to the whole hydrochemical dataset, without separating samples from different locations or supposedly belonging to different water types.This allows that these features, if actually present in the studied system, emerge independently from the statistical analysis.

General hydrochemical features of the Ebro alluvial waters
The general hydrochemistry of the studied portion of the Ebro alluvial-karst aquifer system in the surroundings of Zaragoza city is displayed in Table 2.The observed hydrochemistry for all the points corresponds to Cl-SO 4 -Na-Ca-waters with pH values around 7.2 but varying between 6.5 and 8.6.As a general trend, the highest pH values are measured in surface waters used for irrigation (points IR1, IR2 and IR3 from Table 1), which are also among the points displaying the lowest alkalinity and conductivity values (Table 3).There is a wide Total Dissolved Solids (TDS) variation among different sampling points, roughly ranging between 300 and 6,000mg/l.Such variations are also frequently found among different samples from the same sampling point, as shown by the large standard deviations calculated for some of them (Table 3).

Application of principal component analysis
Principal component analysis has been applied to identify the main variables responsible for the general hydrochemistry of the studied dataset and the possible groups of samples based on their hydrochemistry.After preliminary tests, temperature and dissolved potassium data were excluded from the analysis because their inclusion caused a clear decrease (around 10%) in the percentage of variance explained by the two first principal components without changing substantially the results and interpretations.Therefore, this analysis was finally carried out using the following hydrochemical variables: pH, total dissolved solids and molar dissolved contents of Na, Cl, Ca, Mg, sulfate and alkalinity (as HCO 3-).
The results obtained are shown in Figure 3A and B. Two principal components (PC1 and PC2) were identified explaining about 66 and 20% of the total variance, respectively.Consistently with the Kaiser criterion (Kaiser, 1960), the rest of the components were considered not significant due to their low eigenvalues (<1) and of the low percentage of variance explained by them (below 8%).
In the light of the results, the main processes explaining the variance of hydrochemical data in the Ebro alluvialkarst aquifer waters are related to the increase in the amount of dissolved solids.This is suggested by the association of total dissolved solids and concentrations of Na, Cl, Ca and sulfate as the main variables included in PC1 which, as indicated above, explained almost 70% of the total variance.As shown in Figure 3B, the loadings for all these variables in PC1 are about 0.4.The most probable processes controlling the evolution of these variables in the studied area are the dissolution of halite and gypsum/anhydrite.This issue is further assessed in the following sections.
With regard to PC2, the main variables included are pH and alkalinity (Fig. 3B), with similar absolute loadings in this component of about 0.6.Since the fraction of variance explained by PC2 (around 20%) is much smaller than for PC1, it can be inferred that the influence of processes responsible for alkalinity and pH variations in the system is minor compared to the ones causing salinity increase.Given the reported mineralogy in the alluvial and evaporitic aquifer, the most probable processes responsible for the alkalinity and pH variations in the system seem to be carbonate dissolution (mainly calcite and dolomite) and input/loss of CO 2 (g).The thermodynamic feasibility of such heterogeneous processes in the Ebro alluvial waters will be examined in the next section with the assistance of geochemical calculations.

G e o l o g i c a
(5) Na/Cl (6)  (1) TDS: Total dissolved solids (in mg/l); (2) T: temperature (in o C); (3) Cond: electric conductivity (in µS/cm); (4) Alk: Alkalinity, as HCO 3 -in mg/l; (5) all values in mg/l; (6) molar ratios; (7) dissolved K concentrations were not available for some of the samples taken from the CHE database, and they have been estimated from other samples with similar TDS contents obtained at the same points during the sampling surveys specifically done for this work.

3
Average hydrochemical features for each of the selected sampling points included in this study pond HA1) plot in the region of negative scores for both components.This is consistent with the recharging role of irrigation and with the small degree of interaction with the alluvial and evaporitic materials of the samples included in this group.
ii) Concentrated waters (i.e.waters with total dissolved solids content higher than the mean value for the area, identified as A-type samples; WA, SA and HA sampling points) plot in the region of the positive and largest PC1 values (Fig. 3A).This is consistent with the association of salinity with the PC1.This group seems to include the waters with the greatest degree of interaction with the evaporitic aquifer materials and/or groundwaters and includes, among others, the samples from the three springs (SA2, SA2 and SA3) and some of the samples from one of the sinkhole ponds (HA1) from this study.
iii) Intermediate waters between irrigation and concentrated waters (Fig. 3A).This group includes mostly the relatively dilute samples from B-type wells and boreholes, which plot over the positive region of the PC2 axis and, generally, with values between -2.5 and 1 for the PC1.Their position suggests that these waters are less affected by processes leading to salinity increase than concentrated waters but similarly or even more affected by the processes responsible for the alkalinity variations (PC2 values sometimes larger than for concentrated waters).This group probably represents a variable degree of interaction between the groundwaters of the alluvial and evaporitic aquifers and includes most of the samples from the sinkhole ponds selected for this study.
As explained, the samples corresponding to sinkhole ponds are widely scattered among the three groups described above.In some cases, this dispersion also occurs for different samples from the same sampling point (Fig. 3A).Since sinkhole ponds potentially represent windows in the alluvial aquifer interacting with the underlying evaporitic karst system, their large hydrochemical variability points towards the existence of also variable interaction processes.

Application of hierarchical cluster analysis
For the application of hierarchical cluster analysis to the dataset from the Ebro alluvial-karst aquifer system, two different cluster analyses have been carried out: i) Clustering of samples based on their hydrochemical features, in order to detect possible groupings of sampling points.
ii) Clustering of hydrochemical variables, including dissolved element contents, calculated saturation indexes with respect to the main possible mineral phases reported for the study area (halite, gypsum, calcite, dolomite) and calculated partial pressures of CO 2 (g) in equilibrium with the samples.This will allow identifying the possible controlling geochemical processes in the system (mineral dissolution/precipitation or CO 2 (g) input/loss).
For the clustering of samples, the selected variables were the same ones included in the previous principal component analysis.Hierarchical cluster analysis using Summary of results of Principal Component Analysis (PCA) applied to the most significant hydrochemical variables from the studied system.A) Distribution of sample hydrochemistry in a plot of the first (PC1) vs. the second (PC2) principal components; B) directionality and relative magnitude of influence of the variables used in the PCA on a PC1 vs. PC2 plot, indicated by arrow directions and lengths (respectively).A c t a , 1 1 ( 4 ) , 3 8 9 -4 0 7 ( 2 0 1 3 )  D O I : 1 0 . 1 3 4 4 / 1 0 5 . 0 0 0 0 0 2 0 5 2 Hydrogeochemistry of evaporite Karst in the Ebro valley (NE Spain) 398 Ward's method leads to a dendrogram in which five subbranches are individualized (Fig. 4).As expected from results, of the principal component analysis variables included in the PC1 seem to play a major role in the clustering of those points.Two sample sets are separated at similarities about -44, corresponding to the following main branches indicated in Figure 4: i) Branch 1: dilute groundwaters.From this major branch, two subbranches are separated at similarities about -22: subbranch 1.1, which corresponds to the most dilute waters in the area (most irrigation waters and dilute samples from ponds), and subbranch 1.2, which groups intermediate waters mostly from B-type wells.These subgroups match almost perfectly the corresponding ones in the principal component analysis (Fig. 3A), which provides additional support for the results and of both multivariate analyses.

G e o l o g i c a
ii) Branch 2: concentrated groundwaters, hydrochemically closer to the features of the evaporitic aquifer groundwaters.Two dendrogram subbranches (identified as 2.1 and 2.2 in Fig. 4) are separated at a similarity value around -19, each of them grouping samples with clearly different Cl:SO 4 ratios, as will be explained in the next section with the assistance of ion correlations.Apart from these two subbranches, there is another one (identified as 2.3 in Fig. 4) only grouping the samples from points SA2 and WA7, which are the most clearly related to the evaporitic aquifer in the area east of Zaragoza city.
As described for the principal component analysis results, samples from sinkhole ponds are distributed among three of the five described subbranches (Fig. 4), ranging from concentrated to very dilute waters (even for samples from the same sampling point in the case of sinkhole pond HA1).Their hydrochemistry could be attributed to mixing between highly variable contributions from recharge/irrigation waters interacted with the alluvial materials and groundwater discharge from the underlying evaporitic bedrock.
With regard to the clustering of hydrochemical variables (Fig. 5), the results are in good agreement with the interpretations obtained based on principal component analysis .As displayed in Figure 5, three main branches are obtained in the dendrogram for similarities lower than -21: i) The first dendrogram branch is represented only by temperature (Fig. 5), which is consistent with the poor  Hydrogeochemistry of evaporite Karst in the Ebro valley (NE Spain) 399 correlation obtained between this parameter and the rest of the hydrochemical variables included in this study (Table 4).
ii) The second branch includes most of the measured hydrochemical variables.The saturation indexes with respect to halite and gypsum are included in the same dendrogram subbranches as their corresponding contents of dissolved elements (Na-Cl and Ca-SO 4 , respectively; Fig. 5).This suggests that the dissolution of both minerals, mainly in the evaporitic aquifer according to the mineralogical information, plays a key role on the hydrochemical features of the studied Ebro alluvial-evaporitic system.According to hierarchical cluster analysis results, total dissolved solids contents seem to be more related to gypsum/anhydrite dissolution processes than to halite dissolution, which would be consistent with the predominance of gypsum/ anhydrite in the evaporitic bedrock in this area.With regard to the dissolved contents of Mg, they are linked to the same subbranch grouping total dissolved solids and gypsum/anhydrite dissolution (Fig. 5), which suggests that the dissolution of the main mineral sources for this element is also taking place in the evaporitic aquifer.Potassium dissolved concentrations seem to be weakly linked to the rest of hydrochemical variables in this branch, which is in agreement with their weak correlations (Table 4) and suggests a completely different hydrochemical control of this element in the system.Finally, alkalinity contents and calculated pCO 2 , which are grouped together, are also weakly linked to the rest of variables in this branch.Their association is consistent with the good correlation between these parameters (Pearson correlation coefficient about 0.71) and could be attributed to the key role played by the pCO 2 in the dissolution of carbonate minerals, which must be the main process responsible for the alkalinity contents in the groundwaters.
iii) The third dendrogram branch groups together pH and the saturation indexes calculated for calcite and dolomite.This is not surprising, given the key influence played by the pH in the heterogeneous (i.e.dissolution/ precipitation) and homogeneous (i.e.carbon speciation) processes affecting these mineral phases.
In order to confirm and explain these interpretations and groupings, the processes able to determine the salinity and alkalinity of the Ebro alluvial waters and to determine their hydrogeochemical evolution from recharge waters need to be addressed.This will be explored in the following section with the assistance of ion-ion plots, ion correlations and geochemical calculations.

Ion-ion plots, ion correlations and speciation-solubility calculations
As displayed in Table 4, and as expected from the principal component analysis and hierarchical cluster analysis results, the hydrochemistry of the system is characterized by a good correlation between total dissolved solids values and the dissolved contents of most of the analysed elements (Na, Cl, Ca, Mg and sulfate), yielding for all of them Pearson correlation coefficients above 0.7.
Apart from this, one of the most remarkable hydrogeochemical features that can be deduced for the studied system is the good 1:1 correlation (within ±5% analytical error) between the molar dissolved contents of Na and Cl in most of the studied waters, as shown in Figure 6A.This is also evidenced by the Pearson correlation coefficient between Na and Cl contents, which is around 0.99 (Table 4).Consistently with the hierarchical cluster analysis results and with the mineralogy of the evaporitic bedrock, this observation can be attributed to halite dissolution as the main source for Na and Cl in the system.However, even the water samples with the highest total dissolved solids contents and sampling waters from the evaporitic aquifer are very far from equilibrium with respect to this mineral phase (Fig. 7A).This suggests that the amounts of halite in the sampled portions of evaporitic aquifer are not large enough to attain equilibrium, which is consistent with the predominance of gypsum/anhydrite reported for these materials.

Table 4: Pearson correlation coefficients for hydrochemical variables
SO 4 stoichiometry can be mostly attributed to gypsum and/ or anhydrite dissolution.These phases are among the main minerals described for the evaporitic aquifer.Gypsum is also present in the alluvial deposits, derived from the gypsum bedrock that crops out in the valley margins and from direct precipitation induced by evapotranspiration.In line with this explanation, saturation states with respect to gypsum approach equilibrium with increasing total dissolved solids contents and many of the most concentrated samples from different wells and springs are already equilibrated (within the ±0.22 uncertainty range; Langmuir and Melchior, 1985) with respect to this mineral (Fig. 7B).
The fact that the contents of Na, Cl, Ca and SO4 from most of the samples are linearly located between the values from the recharge waters (irrigation) and from the groundwaters most probably associated with the evaporitic materials (Figs.6A, B) suggests that mixing between these two water types plays a key role in the hydrochemistry of the Ebro alluvial-karst aquifer system.This is perfectly consistent with the groupings of samples inferred from principal component analysis and hierarchical cluster analysis.Furthermore, the observed Na:Cl and Ca:SO 4 correlations and stoichiometries suggest that other possible processes, such as glauberite dissolution or Na/Ca-exchange, generally play a minor role (compared to halite and gypsum dissolution) in the hydrogeochemical behavior of this system.If these processes were of major importance in the whole studied system, they would shift molar Na/Cl and Ca/ SO 4 ratios towards values more clearly different from 1 than generally observed (as displayed, for instance, by the samples from the spring SA3; Fig. 6A, B).Given the large amounts of glauberite reported for the evaporitic succession, the scarce contribution of its dissolution to the groundwater hydrochemistry may be attributed to the preferential circulation of groundwater through areas where it has been replaced by secondary gypsum.Such replacement affects a thick weathering zone as much as 120m deep (Salvany, 2009;Fig. 1B).For the spring SA3, its hydrochemistry is probably affected by glauberite dissolution, since glauberite beds tens of meters thick have been encountered in this specific area at shallow depths (Salvany, 2009).
Good Na:Cl and Ca:SO 4 correlations are also observed in the irrigation waters (Fig. 6A, B).This implies that the same processes of halite and gypsum/anhydrite dissolution are also important in the areas located upstream of the studied area in the Ebro and/or in the Gállego River (Fig. 1A).Consistently, the Ebro River upstream of the studied area traverses outcrops of Paleogene halite-bearing evaporite formation (Falces and Lerín formations), whereas the Gállego River flows through the Zaragoza formation along 40km.
In spite of the already described strong Na:Cl and Ca:SO 4 correlations for the whole dataset, the correlation between Cl and SO 4 dissolved molar contents (Fig. 6C) and, therefore, the Na:SO 4 and Ca:Cl correlations (plots not shown), are much weaker (Pearson correlation coefficients below 0.6; see Table 4).A detailed examination of the plot of Cl vs. sulfate dissolved contents (Fig. 6C) shows two different evolution trends: i) Chloride and sulfate contents increasing steadily and Cl/SO 4 molar ratios between 1.4 and 4.2 (Group 1 in Fig. 6C; see values of the Cl/SO 4 ratio in Table 3).
ii) Low chloride contents (below 20mmol/L) independently of the dissolved sulfate concentrations and Cl/SO 4 molar ratios between 0.3 and 1.2 (Group 2 in Fig. 6C; see values of the Cl/SO 4 ratio in Table 3).
Pearson coefficients for the correlations between the hydrochemical variables included in this study.For the correlations with K, only measured values have been considered The only clear outliers from these groups are the samples from the WA7 point.As can be observed in Figure 7A and B, the two groups differentiated on the basis of their dissolved Cl vs. SO 4 evolution (Fig. 6C) also follow different trends when approaching equilibrium with respect to halite and to gypsum associated with the increase of their total dissolved solids contents.For similar total dissolved solids values, samples from Group 1 are closer to equilibrium with respect to halite but more undersaturated with respect to gypsum than samples from Group 2. Similar results, though with slightly lower saturation indexes, are also obtained with respect to anhydrite.These features suggest that halite and gypsum/ anhydrite dissolution processes have a spatially variable contribution to the hydrochemistry of the Ebro alluvialkarst aquifer system.An examination of the geographical distribution of the samples in both groups (Fig. 1A; Fig. 6C) indicates that samples in Group 1 are all located east of Zaragoza city, whereas samples included in Group 2 correspond to the area west of Zaragoza city (except samples from point WA8).This pattern shows that, as a general trend, the sampled waters in the eastern sector have dissolved larger halite amounts than the ones from the western sector.The most probable explanation for this behavior would be the existence of highest amounts of halite (relative to gypsum/anhydrite) along the flow paths followed by groundwater that arrives at the area east of Zaragoza city.
Apart from gypsum/anhydrite dissolution, the dissolved contents of Ca may also be affected by the dissolution/precipitation of carbonate minerals.The dissolution of these minerals (mainly calcite and dolomite) could account for the trend of increasing alkalinity and Ca contents observed for the samples with less than 0.008mol/l of Ca (Fig. 8A) this could also explain the Ca/SO 4 ratios above 1 observed for the same samples.For the water samples with higher dissolved Ca contents, the roughly constant alkalinity contents observed for any given sampling point and the general trend of decreasing alkalinity for increasing Ca contents for this whole group of samples (Fig. 8A) may be attributed to the attainment of calcite equilibrium and precipitation at low alkalinity values.This can be explained by the presence of increasingly larger calcium concentrations derived from gypsum/anhydrite dissolution (common ion effect).
The apparent oversaturation with respect to calcite (Fig. 7C) and dolomite (Fig. 7D) observed in many Correlation between some of the main hydrochemical parameters in the water samples from the studied system.A) Chloride vs. sodium dissolved concentrations; B) Sulfate vs. calcium dissolved concentrations; C) Chloride vs. sulfate dissolved concentrations, with indication of the two different evolution trends described in the text.A c t a , 1 1 ( 4 ) , 3 8 9 -4 0 7 ( 2 0 1 3 )  D O I : 1 0 . 1 3 4 4 / 1 0 5 . 0 0 0 0 0 2 0 5 2 Hydrogeochemistry of evaporite Karst in the Ebro valley (NE Spain) 402 samples (especially in the ones from sampling points in contact with atmospheric conditions, such as irrigation waters and ponds), may be attributed to CO 2(g) degassing.This is a common situation observed when waters initially equilibrated or even undersaturated with respect to these minerals and with CO 2(g) partial pressures (pCO 2 ) larger than the atmospheric one are put in contact with the atmosphere during sampling or previously (López-Chicano et al., 2001;Auqué et al., 2009 and references therein).Two evidences support this explanation: i) calculated pCO 2 larger than the atmospheric one (around 10 -3.5 atm) in equilibrium with most of the studied waters  (Langmuir and Melchior, 1985), ±0.35 for calcite and ±0.5 for dolomite (Plummer et al., 1990).(Fig. 7E, F), and ii) presence of the largest oversaturation degrees with respect to calcite (Fig. 7E) and dolomite generally associated to samples from surface waters and close to equilibrium with atmospheric pCO 2 .The fact that large pCO 2 are attained in most of the samples points towards the probable CO 2 (g) input to the groundwater system from soils or organic matter degradation, since other possible sources have been neither identified nor reported in earlier works.

G e o l o g i c a
The pH and alkalinity values are also correlated in the system (Fig. 8B), although their correlation is not very strong (Pearson correlation coefficient close to 0.6; Table 4).This can be attributed to the role played by the dissolution and precipitation of carbonate minerals in the control of both variables.
With regard to Mg, the dissolved contents of this element are mainly correlated with Ca and sulfate contents, although this correlation is weaker than the ones described for Ca-SO 4 and Na-Cl (Table 4; Fig. 8C).However, the presence of magnesium sulfates has not been reported for the evaporitic aquifer and its presence in trace amounts in the structure of gypsum/anhydrite (Kushnir, 1980;Lu et al., 1997 and references therein) cannot justify the observed dissolved concentrations in the studied system (Table 3).Therefore, in the light of the reported mineralogy for the evaporitic and alluvial aquifer, the most likely process responsible for the primary release of Mg to the groundwaters is the dissolution of dolomite (CaMg(CO 3 ) 2 ).
Since the highest magnesium contents (between 50 and 100 mg/l) are found in samples associated with the Correlation between some of the main hydrochemical parameters in the water samples from the studied system.A) Alkalinity (as HCO 3-) vs. calcium concentrations, with indication of the two different evolution trends described in the text; B) pH vs. Alkalinity (as HCO 3-); C) Sulfate vs. magnesium dissolved concentrations.A c t a , 1 1 ( 4 ) , 3 8 9 -4 0 7 ( 2 0 1 3 )  D O I : 1 0 . 1 3 4 4 / 1 0 5 . 0 0 0 0 0 2 0 5 2 Hydrogeochemistry of evaporite Karst in the Ebro valley (NE Spain) 404 evaporitic aquifer (A-type wells and springs; Fig. 8C; Table 3), this process seems to have taken place more prominently in these materials.This would explain the observed correlation of Mg with Ca and sulfate and would be consistent with the hierarchical cluster analysis results.Two complementary processes may account for the large dolomite dissolution in these A-type groundwaters: i) dedolomitisation triggered by gypsum/ anhydrite dissolution (e.g.Plummer et al., 1990), which is supported by the fact that these waters are close to calcite and dolomite equilibrium (Fig. 7C, D, E) independently of their closeness to gypsum equilibrium (Fig. 7B), and ii) presence of high pCO 2 values (above 10 -2 atm; Fig. 7E, F) in equilibrium with them, which is probably caused by CO 2 input to recharge waters from biological activity in soils.

G e o l o g i c a
On the contrary, samples displaying a hydrochemistry apparently more linked to the alluvial aquifer (e.g.most samples from B-type wells), generally show lower Mg concentrations (below 50mg/L; Table 3).This may be due to the existence of lower amounts of dolomite in the alluvial aquifer and/or to the lower influence of dedolomitisation.
For the dissolved molar K concentrations no significant correlation with Cl, SO 4 or other dissolved molar contents is observed (Table 4), which is consistent with the weak linkage of the contents for this element to other ions in the hierarchical cluster analysis.This behavior is not surprising, since the main sources for potassium in these groundwaters are probably aluminosilicate phases or even agricultural fertilizers, whose distribution and dissolution behavior must be very different from halite, gypsum/ anhydrite and even carbonate minerals.Moreover, there are not significant differences in K dissolved concentrations between different water types, which prevents from knowing if the processes responsible for K release are preferentially linked to the alluvial aquifer or the evaporitic materials.

CONCLUSIONS
The application of multivariate analyses, ion correlations and speciation-solubility calculations to a comprehensive hydrochemical dataset from the Ebro alluvial aquifer and the underlying karstified salt-bearing evaporitic bedrock has allowed the identification of some of the main features and processes of the system.The results obtained are the first step to shed light on the inner-workings of the alluvial-karst aquifer system and to improve our understanding of the karstification processes behind sinkhole development in the area.
As suggested by principal component analysis, hierarchical cluster analysis, ion correlations and geochemical calculations, the hydrogeochemistry of the system seems to be mainly determined by the variable contribution of gypsum/anhydrite and halite dissolution, showing a wide spectrum whose end members correspond to relatively fresh recharge waters (i.e.waters used for irrigation) and highly concentrated groundwaters related to the discharge of water from the evaporite karst aquifer in the Ebro alluvial aquifer.In any case, whereas equilibrium with respect to gypsum is attained by some of the groundwaters sampled in the evaporitic aquifer, halite equilibrium is far from being reached even for the most concentrated groundwaters.This may be attributed to the existence of low amounts of halite (relative to gypsum/ anhydrite) interacting with the sampled groundwaters in the evaporitic bedrock.Moreover, the relative contribution of gypsum/anhydrite and halite dissolution to the observed hydrochemistry is clearly spatially variable.As a general trend, the hydrochemistry of groundwaters from the Ebro alluvial aquifer and from the underlying evaporitic materials west from Zaragoza city seems to indicate lower interaction with halite than the ones sampled east of Zaragoza city.Another interesting mineralogical and geochemical observation is related to the minor role played by glauberite dissolution on the observed hydrochemistry, even though this mineral represents a significant proportion of the unweathered bedrock.This can be due to the extensive replacement of glauberite by secondary gypsum in most of the evaporitic deposits interacting with the groundwaters.
Apart from these karstification processes, the hydrochemistry is also influenced by the dissolution of carbonate minerals (dolomite and, possibly, calcite) triggered by the input of CO 2(g) , which is probably produced by pedogenic processes in this area with ubiquitous agricultural activities.These processes seem to be responsible for the evolution of alkalinity and pH buffering.Moreover, they may also increase dolomite dissolution in the evaporitic aquifer, which is suggested by the larger Mg concentrations observed in the groundwaters linked to the karst aquifer compared to the recharging irrigation waters.Even though dolomite dissolution may also occur in the alluvial aquifer, this process seems to be more relevant in the evaporitic materials and to be also related to dedolomitisation triggered by gypsum/anhydrite dissolution.In the case of calcite, both dissolution and precipitation of this mineral phase are probably taking place in different portions of the studied system, although this issue should be addressed further in future studies.
The hydrochemistry observed at the sampled perennial sinkhole ponds displays a large variability, from weakly evolved waters very similar to the irrigation waters, to others with geochemical features closer to the ones observed in deep wells abstracting water from the evaporite karst aquifer, which seem to be associated with discharge areas of the karst aquifer in the valley bottom alluvium.This spatial relationship may be useful to assess sinkhole susceptibility in the area.Although the development of sinkholes in the discharge areas of carbonate karst aquifers is quite unusual (Salvati and Sasowski, 2002;Tuccimei et al., 2005), it seems to be more likely in karst systems developed in halite-bearing evaporites where discharging groundwaters are not saturated with respect to halite (e.g.Land, 2003;Yechieli et al., 2006)., 1 1 ( 4 ) , 3 8 9 -4 0 7 ( 2 0 1 3 )  D O I : 1 0 . 1 3 4 4 / 1 0 5 . 0 0 0 0 0 2 0 5 2 Hydrogeochemistry of evaporite Karst in the Ebro valley (NE Spain) FIGURE 1 Images of some of the sampled permanent sinkhole ponds.A) Ojo del Fraile, a shallow lake (around 240m of maximum length) developed in a sagging sinkhole with nested collapse dolines (sampling point HB1).B) Torre del Chocolatero collapse sinkhole and pond (sampling point HB2).A maximum water depth of 6.6m was measured in April 2009 in this pond (around 55m of maximum length).C, D) Ojos de Matamala flooded collapse sinkholes.This roughly circular pond (sampling point HA1) is 30m in diameter and 3.5m deep (as measured in April 2009).

FIGURE 2
FIGURE 2 FIGURE 3 Simplified results from Hierarchical Cluster Analysis (HCA) applied to samples and sampling points of the studied Ebro alluvial waters.For a better visualisation of results, only clustering for similarities above -15 and sampling points (instead of individual samples) are shown.Sampling point names surrounded by a rectangle are the ones for which every single sample falls into the displayed grouping, whereas sampling point names without a rectangle indicate than some of the samples from those points belong to different cluster branches.Branches and subbranches have been numbered and are described in detail in the text.
FIGURE 4 FIGURE 5 FIGURE 6 Correlation between calculated saturation indexes with hydrochemical parameters in the water samples from the studied system.A) Halite saturation index vs.Total Dissolved Solids (TDS), with indication of the two different evolution trends described in the text; B) Gypsum saturation index vs.TDS, with indication of the two different evolution trends described in the text; C) Calcite saturation index vs.TDS; D) Dolomite saturation index vs.Calcite saturation index; E) log[pCO 2 (g)] vs Calcite saturation index; F) log[pCO 2 (g)] vs. TDS.The shaded areas indicate the ranges of uncertainty for the calculated saturation indexes: ±0.22 for gypsum

FIGURE 7
FIGURE 7 FIGURE 8 l o g i c a A c t a

TABLE 1 Table 1 : General features of sampling points (modifiable) 1
(1) hydrochemical data obtained from the Ebro Waters Authority (CHE); (2) hydrochemical data obtained from sampling carried out for this study; (3) according to data from CHE water points database (IPA) and/or provided by drilling and engineering geology companies; (4) according to waterlevel data at wells and to alluvial aquifer thickness; (5) elevation data obtained from the CHE Digital Elevation Model (DEM); n.a.: information not available; -: information not applicable for this type of sampling point.P. A c e r o e t a l .