A geometrical method for quantifying endmembers' fractions in three‐component groundwater mixing

We present a new geometrical method capable of quantifying and illustrating the outcomes of a three‐component mixing dynamics. In a three‐component mixing scenario, classical algebraic equations and endmember mixing analysis (EMMA) can be used to quantify the contributions from each fraction. Three‐component mixing of natural waters, either in an element–element plot or by using the EMMA mixing subspace is described by a triangular shaped distribution of sample points where each endmember is placed on an apex, while each side corresponds to the mixing function of the two endmembers placed at the apex, considering the third endmembers' contribution equal to zero. Along each side, the theoretical mixing fractions can be computed using mass balance equations. Samples with contributions from three endmembers will plot inside the triangle, while the homogeneous barycentric coordinate projections can be projected onto the three sides. The geochemistry observed in the mineralized Ferrarelle aquifer system (southern Italy) results from three‐component mixing of groundwater, each with diagnostic geochemical compositions. The defined boundary conditions allow us to parameterize and validate the procedures for modelling mixing, including selection of suitable geochemical tracers.


| INTRODUCTION
Mixing between different sources of water is one of the most common and important hydrological process in nature (e.g., Barros Grace et al., 2008;Domenico & Schwartz, 1997;Ram on et al., 2021;Stegen et al., 2016;Wilson et al., 2016). The definition of the endmembers of the mixing allows to understand the hydrodynamics and the chemical processes.
The results of these mixing processes can sometimes influence water quality limiting water access due to the presence of lower quality water mixtures. The worsening of the original groundwater composition could be related to natural contaminations (uprising of deep geothermal fluids rich in arsenic (As), iron (Fe), or other trace metals or hydrocarbons, Bundschuh & Maity, 2015;Donato et al., 2021;Tassi et al., 2009) or due to leaching of shallow aquifers contaminated by human activities (pesticides and fertilizers in agricultural areas or wastewaters in urban areas, Li et al., 2021;Mateo-Sagasta et al., 2017). As a result of these processes, the final composition of a mixed parcel of water can be complex and depends on the relative concentrations and percentage of each mixing component.
The investigation of hydrological processes related to the groundwater mixing dynamics can be especially important in complex hydrogeological systems, where vertical and/or lateral groundwater exchanges occur between different aquifers, or where the lack of data limits the understanding of deep aquifer dynamics.
The effect on solute concentrations during mixing are traditionally described using mathematical models (e.g., Christophersen & Hooper, 1992;Mazor et al., 1993;Schramke et al., 1996). Using these methods, algebraic modelling can be used to quantitatively determine the fractions of original water (endmember) and added components that were contributed to the sample mixture. More recently developed mixing models (Hooper, 2003;Liu et al., 2017) are based on computational linear algebra, known as end member mixing analysis (EMMA). The EMMA method is widely employed for exploration of complex hydrological systems (James & Roulet, 2006;Katsuyama et al., 2001;Kronholm & Capel, 2014;Liu et al., 2008), as well as groundwater systems in which chemical reactions occurring in a diverse range of hydrogeological settings are of great significance (Pelizardi et al., 2017).
Robust detection and modelling of groundwater mixing processes will allow the improvement of hydrogeological conceptual modelling and improve the understanding of groundwater geochemistry, which are important first steps in more elaborate investigations on water (e.g., water budget calculation, evaluation of potential contamination from external sources). The ability to improve conceptual models is an important tool for developing proper management strategies, especially in mineralized aquifers tapped for bottling activity, where the healthiness and the constancy of the groundwater chemical composition must be guaranteed (Europe Directive 2009/54/EC). When endmembers mix conservatively the mixture, geometrically, will lie in the convex set where the vertices are the three components (Christophersen & Hooper, 1992 and therein references). In case of a three component mixing, the set corresponds to a triangle on the plane defined by the three endmembers. This triangle can act as reference triangle in a homogeneous barycentric coordinate system.
In the present investigation we test a new and fast geometrical method to quantify the relative mixing fractions of three endmembers directly from the bidimensional mixing space. This method can be applied both for the endmember fractions estimations based on the simple mass balance equations and for the multivariate EMMA method.
The proposed method was tested on the Ferrarelle groundwater system (FGS), a complex mineral system strictly governed by groundwater mixing (Cuoco et al., 2020). Despite the complex and peculiar hydrogeological setting, the studied system is supported by an extensive bibliographic knowledge (Sacchi et al., 2021 and reference therein) and a robust monthly chemical dataset. Specifically, the water chemistry reflects a balance between recently entrained volcanically influenced waters and carbonate aquifers (Cuoco et al., 2010). Upflow of CaCO 3 and CO 2 rich groundwater occurs from the deep confined carbonate aquifer through the fault systems, which greatly influence the sedimentary bedrock and allow mixing with less mineralized groundwater in the volcanic aquifer (Corniello et al., 2015;Cuoco et al., 2017;. Cuoco et al. (2020) and  argued for the presence of a third groundwater component in the FGS, which was hypothesized to result from mixing by lateral inflow from the Mt. Maggiore region. This third component is dominantly alkali earth-bicarbonate-type water with low total dissolved solids. Thus, Cuoco et al. (2020) (Boncio et al., 2016;Cosentino et al., 2006;Giordano et al., 1995). The Meso-Cenozoic carbonate units crop out both in the Mt. Maggiore (D'Argenio & Pescatore, 1962) and in the other mountain ridges that surround the Riardo Plain. Following deformation during the Apennine orogeny, the carbonate deposition ended and gave way to the deposition of marl and synorogenic flysch units. The flysch units are not deposited or preserved continuously as both processes were controlled by the presence of localized graben structures throughout the Riardo Plain (Saroli et al., 2014), whereas horsts are zones of either non-deposition or erosion.
The volcanic sequences within the Riardo Plain can be summarized as a sequence of loose layered pyroclastic and reworked volcanoclastic deposits from the Roccamonfina volcano (Viaroli, Lotti, et al., 2019) with alternating deposition of the Brown Leucitic Tuff F I G U R E 1 (a) Scheme of the regional aquifers in the study area (modified from Boni et al., 1986;Celico, 1978;Lancia et al., 2020); (b) Location of the available boreholes and the monitored wells in the Ferrarelle mineral water claim area; (c) Hydrogeological cross sections of the study area. The ID codes of sampled wells (in blue) are taken from Cuoco et al. (2020) unit, which is a lithic tuff composed of an ash matrix (Luhr & Giannetti, 1987) (Figure 1c). Petrographic analysis carried out on cutting samples from boreholes in the Ferrarelle area (Università degli Studi Roma Tre, 1996, technical report) indicate the presence of several clay minerals (i.e., Halloysite 10 Å, Smectite, Chabazite and Illite) derived from the weathering of the volcanic glasses.
The basal aquifer of the Roccamonfina volcano shows a radial drainage pattern towards the plains that surround the volcano edifice.
Groundwater in this aquifer is mainly discharged through the Savone gaining stream with a mean discharge of $0.5 m 3 /s (Boni et al., 1986).
The Mt. Maggiore aquifer feeds the Triflisco Spring, located at the southern side of the carbonate structure at an elevation of $30 m a.s. l. with a mean discharge of $2.8 m 3 /s (Boni et al., 1986) (Figure 1a).
In the Riardo Plain, the hydrogeological framework is quite complex and heterogeneous. In the western Riardo Plain, the absence of intermediate aquiclude units allows for the upwelling of deeply source mineralized fluids along faults affecting the sedimentary basement, which leads to mixing between the shallow volcanically influenced groundwater (N component) and the more deeply sourced and highly mineralized carbonate aquifer (D component), the latter of which is responsible for creating the Ferrarelle mineral springs (Cuoco et al., 2010(Cuoco et al., , 2020Sacchi et al., 2021). The mineralization of groundwater varies locally (well to well) with groundwater chemistry varying according to the proximity of faults. While local chemistry is governed by discrete localized zones of fault-controlled fluid migration, the vast majority of data from the study area fall along an apparent twocomponent mixing line between volcanically influenced groundwater and the more deeply sourced and highly mineralized carbonate aquifer (Cuoco et al., 2020). Recent publications, Cuoco et al. (2020) identified the presence of an additional contribution in the mineral water mixing, corresponding to a lateral inflow from the phreatic carbonate aquifer of Mt. Maggiore (F component), however the relevance of this contribution in the groundwater mixing is not yet defined.

| Field measurements
Monitoring and sampling activities were performed monthly between October 2017 and January 2019 on eight monitoring wells placed within Ferrarelle™ company's water concession in the Riardo Plain (blue dots in Figure 1b). Groundwater depth was measured in the piezometers where the static condition of the aquifer was well established. Groundwater depth was also measured daily in two monitoring wells (N and F in Figure 1b) using automatic probes (Schlumberger DIVER) during the same period ($18 months). The N and F monitoring wells tap the volcanic aquifer in the Riardo Plain and the carbonate non-mineralized aquifer of Mt. Maggiore, respectively, in order to evaluate the relationship between the aquifers during long-term monitoring. On average, the hydraulic heads measured in the mixing area were $100-115 m a.s.l. (Figure 1c).
Field measurements were performed following the procedures described in Cuoco et al. (2020). Temperature, pH and electric conductivity (EC, with automatic compensation to 20 C) were measured in each surveyed point using portable probes (WTW pH/EC 340i).
Precision and accuracy were determined for pH and EC by measuring relevant WTW certified standard solutions, the results were better than 1%. Redox potential (ORP) was measured in the field using a Hanna HI991002 meter with precision better than 5%.

| Chemical analysis
Major elements were analysed by ion chromatography on unacidified (F,Cl,Br,NO 3 and SO 4 ) and acidified (Na, K, Mg and Ca) samples. Bromide resulted always below detection limit. The total alkalinity (as HCO 3 ) was determined in the field (Kartell™ technotrate digital burette) by titration using standardized 0.1 M HCl (Merck™) against methyl orange indicator. All samples were filtered in the field using 0.45 μm Minisart sterile cellulose acetate membrane filters and then an aliquot was acidified ($1%) with ultra-pure Merck HNO 3 . Precision and accuracy for the IC analyses was tested by measuring certificated reference materials, confirming analytical error less than 7%. All reported water analyses show charge balances of less than ±5%.
Trace metal analysis were carried out using Agilent 7500ce-ORS ICP-MS with collision cell technology according to the procedure described by Cuoco, Darrah, Buono, Verrengia, et al. (2015). Internal standards (constant concentration of 89 Y and 159 Tb) were used to monitor the instrumentation drift. External precision was monitored using an internal spike standard, which indicates errors of less than 10% for all reported analytes. All the trace metal analytes reported herein exceeded the limits of detection and quantification as calculated according to Long and Winefordner (1983) by at least an order of magnitude.

| Theoretical elements
In natural water systems, two or more sources of water can mix and change the resulting chemistry. Mixing can be described by an array of mathematical relationships (e.g., Albarede, 2009;Rutherford, 1994), including algebraic mass balances based on selected chemical variables (i.e., major ions, key elements). During mass balance approaches, the set of variables and equations are defined and used to determine mixing models that are used to quantitatively describe the hydrodynamics of groundwater systems.
The mixing between two different parcels of water is defined as where A is the concentration of the conservative element, which is now defined as a 'tracer', A mix , A α and A β are respectively the concentrations of A in the mixed water, in α and β endmembers. If B is a second, different tracer, the mass balance can be also written: Considering that the sum of the fractions must be = 1: substituting Equation (3) in (1) and (2), and combining (1) and (2), the following relationship is obtained: Relation (4) (1) and (2) and plotted on the line of the mixing function.
In a 'three-component mixing' scenario using the mass balance approach, three different water volumes, named α, β and γ, are assumed to mix each together. The mass balance equations are then defined by three fractions: The mass balance needs to include two tracers (A mix and B mix ) and can be solved as follows: The algebraic solutions of (5), (6) and (7) give numerical values of the three fractions in the resulting water (Genereux, 1998;Liu et al., 2017). The geometric representation of the three-component mixing in a binary plot corresponds to a triangle (i.e., ternary diagram) ( Figure 2b) whose vertices are defined by the coupled values of A and B tracers detected in the three endmembers. On the triangle surface, the water samples which are produced by the mix of the three endmembers can be plotted (Renner, 1988).
Fixing one of the three components equal to 0, the problem becomes a two-component mixing calculation ( Figure 2a). As consequence, each side of the triangle can report the theoretical fractions of mixing related to the adjacent components ( Figure 2b). From Equation (5) it is given by: Considering Equation (8), the concentration C of a third conservative tracer in a generic solution b produced by the mixing of α and β endmembers is given by the following mass balance equation: F I G U R E 2 (a) A two-component mixing on an element-element plot. The theoretical fractions can be fixed and the corresponding coordinates can be computed through mass balance equations with the goal of reporting the references on the mixing line. (b) The plot shows a three-component mixing triangle where each endmember is placed on a different vertex. The theoretical mixing fractions has been computed by fixing the third endmember as equal to zero. A generic point P produced by a three-component mixing (P A and P B are given by chemical analyses) will be detected in a plot by the homogeneous barycentric coordinates systems, which allows one to show the quantification of each endmember on the corresponding side similarly, from Equations (9) and (10) respectively: In tracer-tracer plots, these equations are expressed in vectorial form (bold terms of the eq. stand for 'vector'), with each concentration expressed by two-components related to tracers A and B.
Each side of the ternary diagram shown in Figure 2b reports the theoretical fractions of mixing computed considering the third endmember is equal to 0. For example, a water solution P, produced by the three-component mixing of α, β and γ, must fall within the ternary diagram ( Figure 2b) and the fractions of each endmember are given by the homogeneous barycentric coordinates related to the theoretical fractions reported on each side of the ternary diagram: being f α þ f β þ f γ ¼ 1 (Lidberg, 2011). The convex combination of the theoretical fractions on the ternary diagram detects the relative concentration of the α, β and γ endmembers which form the P solution.
In the example of the Figure 2b the position of P is defined through coordinates given by the triple of numbers (f α : f β : f γ ) in the specific example is (0.3: 0.2: 0.5) in a given triangle αβγ.
If we imagine the f α , f β and f γ as masses, β and γ will balance at the point X ¼

| Determination of the suitable tracers
The whole dataset of chemical analysis is reported in the dataset stored in the specific repository , while a data summary for all monitored wells is shown in Table 1. In the dataset, each monthly sampling point has been identified with a lower-case letter of Latin alphabet.
The suitability of each chemical variable is evaluated to determine its effectiveness as a potential tracer. The chemical variables involved in the fraction computations for each tracer must respect important constraints in order to get a reliable mixing model (Elsenbeer et al., 1995;Inamdar, 2011;Liu et al., 2017) These constraints include: 1. The tracers in the mixing process must exhibit a linear relationship; 2. The tracers must behave conservatively in aqueous solutions; 3. The components have significantly different concentrations for at least one tracer; 4. The concentrations of each tracer in each endmember must be consistent throughout the investigated area or over the time scale considered for the mixing model.
The mixing process related to the dilution of the more concentrated endmember produces linear correlations among water solutes (Hooper, 2003). In the correlation matrix, reported in the dataset   Note: D, N and F are the endmembers, the other groundwater wells have the same ID as presented in Cuoco et al. (2020). Abbreviations: nd, not detected; RSD, relative SD. , most of the chemical variables are well correlated, indicative of mixing, as demonstrated by Cuoco et al. (2020). This statistical relationship is shown in Figure 4a,b where the EC versus HCO 3 (Figure 4a) and Na + versus K + (Figure 4b) display highly significant correlations (R 2 = 0.99 and 0.95 respectively, p < 0.001). Similar linear correlations (R ≈ 0.90) have been found among the most of highly mobile elements: alkali elements with exception of Rb and Cs, alkali earth elements, HCO 3 , B and EC.
In order to behave conservatively, a variable must not be affected by any chemical reaction. In the investigated system the chemical reactions which can occur are: (i) the re-equilibration of water solution with the CO 2 degassing; (ii) redox re-equilibration due to the mixing between the low ORP deep water and the more oxidant one in the shallow aquifer; (iii) hydrolysis of minerals/volcanic glass and (iv) ion adsorption on the host rock and on chemical precipitates.
According to calculations performed using PHREEQC (Parkhurst & Appelo, 1999), waters with EC >2000 μS/cm are saturated (saturation index, SI > 0) with respect calcite and those with EC >2500 μS/cm are saturated with respect dolomite (see Table 1).  (Drever, 1997;Kehew, 2001;Langmuir, 1997). As a free gas, carbon dioxide migrates in the water column differently than when dissolved in solution and therefore, the chemical equilibrium is not exclusively related to mixing processes of distinct water parcels.
Because of this lack of conservative behaviour, Ca, Mg, HCO 3 and pH (i.e., [H + ]) are thus not considered suitable tracers. EC is strongly correlated to these solutes , and thus is also dismissed as unreliable because of its non-conservative behaviour.
The dissolved oxygen concentration is the main geochemical variable that influences ORP (Stumm & Morgan, 1996). When oxygen levels increase after mixing, Fe and Mn-rich form nearly insoluble oxyhydroxides and precipitate on the surfaces of aquifer minerals (Ahmad et al., 2019); this process is observed in D samples. Here, we observe an inverse correlation (R = À0.83) between Fe and Mn relative to ORP  indicating that both species are likely lost to precipitation after oxidation, as would be expected. Some tracers, including Ba, Sr and As (R $ 0.9), also correlate strongly with Mn and Fe. These trends suggest that these elements are also likely absorbed onto or co-precipitate with Fe/Mn precipitates (Dowling et al., 2002;Rashad et al., 2008) as we previously suggested in this region (Cuoco et al., 2020) and as it has recently investigated in depth . Ion exchange processes may also affect the chemical composition of water solutions (Appelo & Postma, 2005). When groundwater flows into volcanic successions, it is exposed to exchange sites present on the rocks and on Fe and Mn precipitates. The (Na + K)/(Mg + Ca) is commonly used to determine if the adsorption-desorption reactions have altered the chemical composition of the major ions. If the exchange reactions did not affect the composition of dissolved F I G U R E 4 Linear relationship between HCO 3 versus EC (a), and K versus Na (b) as consequence groundwater mixing among endmembers with different in ionic concentrations and CO 2 abundances ions, the ratio in the mixing dynamics would be constant, and therefore the metals concentrations exclusively depending on the dilution effect.
In the D samples, the mill-equivalents ratio (Na + K)/(Mg + Ca) was approximately equal to 0.1-0.2 (the linear correlation of Na + K versus Ca + Mg displayed R 2 = 0.78, suggesting a quite constant angular coefficient, i.e., ratio value) and remained nearly constant throughout the sample set with the exception of samples from unmixed N endmember (0.9-1.0), which reflects the specific ratio of the shallow volcanic aquifer. Therefore, the major alkali ions Na and K are apparently not strongly affected by ion exchange dynamics and can be considered effective as tracers for mixing fraction quantification.
Conversely, trace alkali and alkali earth metals were not considered reliable tracers because their concentrations are 1000 times less concentrated than major ions and thus the effect of interactions with exchange sites may more readily affect their concentration in solution.
Other trace elements such as Li, B, Ba, Cs were not constant in the D aquifer during the monitoring period. The relative SD (RSD) of major elements (Table 1) were all less than 5%, while the trace elements differed by between 10% and 50%. Thus, we would not consider trace elements, with the exception of V, when constructing mixing models.
The last step in the selection of suitable tracers is an evaluation of significant differences in the endmember's concentrations. This decision-making process can be led by observing the statistical distribution of the chemical variables in the dataset. If the concentration of a tracer is similar in all the endmembers, the solution produced by the mixing will not have a significantly different concentration. As a result, tracer behaviour will mimic a random variable and likely display a similar statistical distribution to a normal frequency function. Figure 5 shows the histograms of Na, K, V and Cl concentrations. Note that Cl, which is widely considered as one of the best tracers for mass balances in water mixing due to its highly conservative behaviour in solution (Harkness et al., 2017;Hendry et al., 2000;Warner et al., 2012), follows a normal distribution (Shapiro-Wilk test, p = 0.071); thus, we conclude that Cl is not a suitable tracer in this current study area. Similarly, temperature, pH and F are also normally distributed and not suitable tracers (p > 0.054).

| Characterization of endmembers and related tracers
We consider three endmembers, based on data from the D, N and The N endmember corresponds to the groundwater hosted in the volcanic deposits and it is characterized by low mineralization ($430 μS/cm), the highest ORP (110 mV) observed in this study and low CO 2 concentrations (0.004 mol/kg). We anticipate that the mixing process occurred where more deeply sourced fluids migrated along the fault planes as depicted by the linear relationship (Figure 4) between the mineralization level (EC) and HCO 3 (and indirectly also the CO 2 levels in solution) (Drever, 1997). The linear function is due to mixing, but when one only considers the variables in the binary plots of Figure 4, the third endmember (F) is not distinguished from the N endmember.
The The D endmember is dominated by Ca HCO 3 , but shows the highest concentrations of major ions, alkali and alkaline earth elements and trace elements due to water-rock interactions with the Roccamonfina volcano deposits (Cuoco et al., 2010(Cuoco et al., , 2020Viaroli et al., 2018). N is a Na K HCO 3 type water, typical of groundwater in volcanic aquifers from central Italy (Gambardella et al., 2005). All

| Fractions quantification and analytical representation
After the selection of K, Na and V as tracers, the mixing fractions have been computed using the algebraic solution of the mass balance equations (Liu et al., 2017;Ogunkoya & Jenkins, 1993). For a threecomponent mixture, two tracers are necessary in the calculation. The uncertainties have been computed by Gaussian error propagation (Genereux, 1998). The full results are reported in stored dataset  and summarized the Table 2. The mixing fractions have been obtained both by cross plots of Na versus V and K versus V, whereas the cross plots of Na versus K have not been considered due to the similarity in concentrations between the endmembers (Table 1).

The variability of the calculated fractions was evaluated through
Mann-Whitney rank sum tests (α = 0.05) and were not statistically significant different (p > 0.091) suggesting that the quantifications obtained through Na-V and K-V couples gave similar mixing estimates. In the Figure 6a  The EMMA is a mathematical method based on multivariate statistics, used to identify and quantify the endmembers contribution in complex hydrological systems Christophersen & Hooper, 1992;Hooper et al., 1990). EMMA methods are applied to the whole dataset without any previous supposition about the endmembers. It is only necessary to choose the effective tracers, which will subsequently be tested as part of the procedure (Hooper, 2003). Here, we tested the EMMA method using the same three geochemical tracers (i.e., Na, K, V).
The principal components have been extracted after the data normalization (Joreskog et al., 1976). can be defined as such: where U1 and U2 are the two vectors defining the mixing subspace U, U1 p and U2 p are the coordinates of p (generic sample point whose chemical composition is the produced by the three-component mixing) and D, N and F subscripts detects the coordinates of the three endmembers. Combining Equations (15)- (17), the mixing fractions can be obtained as follows: Mixing fractions computed through the EMMA method have been reported in Cuoco, Viaroli, Darrah, et al. (2021) and summarized in the   Another issue which cannot be neglected for a correct mixing model assessment, is the statistical distribution in the dataset of the tracers. Szulczewski and Jakubowski (2018) observed that a dataset produced by mixing follows a multimodal distribution, as it is confirmed in the Figure 5 for the Na, K and V histograms. Cl, which is usually assessed as a powerful tracers due to its conservative behaviour, in the study system ti shows a normal distribution and consequently the fractions computed using Cl, both using the simple algebraic and  Detailed differences between each monitoring point may also depend on the structure of the groundwater well (e.g., location of the screens, the mean pumping rate, etc.) and on the local structural settings. The available hydrogeological and geochemical conceptual models, coupled with the large dataset of the Ferrarele area, promote this site as an optimal experimental site for scientific research on mixing hydrological processes involving mineral waters.

| CONCLUSIONS
A new method to quantify three endmember mixing was successfully applied in a mineralized groundwater system in southern Italy. The The results improve also the knowledge on the contribution of shallow aquifers, which should be necessary taken in account in a sustainable exploitation plan.

DATA AVAILABILITY STATEMENT
Datasets for this research are available in these in-text data citation references: Cuoco, Viaroli, Darrah, et al. (2021)