A multicomponent geothermometer for high-temperature basalt settings

For successful geothermal reservoir exploration, accurate temperature estimation is essential. Since reservoir temperature estimation frequently involves high uncertainties when using conventional solute geothermometers, a new statistical approach is proposed. The focus of this study is on the development of a new multicomponent geothermometer tool which requires a significantly reduced data set compared to existing approaches. The method is validated against reservoir temperature measurements in the Krafla and the Reykjanes geothermal systems. A site-specific basaltic mineral set was selected as the basis to compute the equilibrium temperatures. These high-enthalpy geothermal reservoirs are located in the neo-volcanic zone of Iceland where the fluid temperatures are known to reach up to 350 °C at a depth of 2000 m. During ascent, the fluid composition is prone to changes as well as possible phase segregation due to depressurization and boiling. Therefore, to reduce the uncertainty of temperature estimations, reservoir temperature conditions are numerically reconstructed with sensitivity analyses considering pH, aluminium concentration, and steam loss. The evaluation of the geochemical data and the sensitivity analyses were calculated via a numerical in-house tool called MulT_predict. In all cases, the temperature estimations match with the in situ temperatures measured at Krafla and Reykjanes. The development of this method tends to be a promising and precise tool for reservoir temperature estimation. The developed methodology is a fast and easy-to-handle exploration tool that can be applied to standard geochemical data without the need for a sophisticated gas analysis yet obtaining very accurate results.

. These approaches use the temperature dependence of the saturation of mineral phases (e.g. silica) or certain cation ratios in the fluids. The measured concentrations of these fluid constituents are then directly linked to a reservoir temperature. The fundamental assumption of geothermometry is the overall chemical equilibrium of the fluid and the reservoir rock. Secondary processes may change the fluid composition and hence, the equilibrium while migrating to the earth surface. These variations can result in large uncertainties for the reservoir temperature determination using solute geothermometers (Nitschke et al. 2018). The more recently developed multicomponent geothermometry evaluates the equilibria of multiple mineral phases . Numerical geochemical speciation codes facilitate this evaluation based on a large number of minerals, which leads to a statistically more robust method. Spycher et al. (2014) proposed a pre-selection of minerals representing the site-specific reservoir rocks to enhance accuracy. Corrections were established to overcome interferences from secondary processes such as dilution, boiling, and mixing of fluids affecting the temperature estimation (Cooper et al. 2013;Peiffer et al. 2014;Spycher et al. 2014). These methods need an additional gas analysis for precise temperature estimations. Thus, Nitschke et al. (2017) introduced a method to reconstruct in situ conditions of the reservoir temperature by varying sensitive parameters, especially pH and aluminium concentration as well as steam loss, to further reduce the uncertainty of equilibrium temperatures.
The goals of this study are to refine and to validate the existing specific multicomponent approach according to Nitschke et al. (2017) and to expand it towards a highprecision exploration tool. This study devises a basalt-specific mineral set including secondary mineral phases for global application to basaltic stratigraphy. For the validation, geochemical data and in situ temperature measurements from basalt-hosted geothermal systems, Krafla and Reykjanes, are used. These are high-enthalpy systems with near-boiling reservoir fluids and, thus, the effect of steam loss between the reservoir and the liquid sampled at the well-head has to be considered. Krafla hosts dilute meteoric fluids , and Reykjanes a more saline reservoir fluid partially originating from seawater (Arnórsson 1978). To validate the method, the temperature estimations are compared with direct in situ temperature measurements of the wells published by Guðmundsson and Arnórsson (2002) for Krafla, and Óskarsson et al. (2015) for Reykjanes. The advantage of the validated method is the frugality in terms of input data. High-accuracy temperature estimations can be achieved based on standard fluid analyses and do not require comprehensive high-end fluid and gas analyses, which are required for other solute multicomponent geothermometer approaches.

Method and data
The basis of the method is a standard fluid analysis comprising major cations and anions as well as aluminium concentration and pH. The water analysis is used to calculate equilibrium conditions between the dissolved constituents in the geothermal fluid and the reservoir rock minerals. For identification of the reservoir temperature conditions, sensitive parameters have to be evaluated statistically. The application relies on the following general assumptions: (i) the reservoir and the geothermal fluid are in equilibrium. Therefore, the ion activity product of a mineral phase equals its thermodynamic equilibrium constant. (ii) A temperature-dependent reaction between the host rock and the water leads to a specific amount of dissolved solids in the fluid phase.
The equilibrium reaction is based on the law of mass action. The state of the dynamic equilibrium between the reactants is expressed in terms of the saturation index SI (1) with IAP being the ion activity product and K being the temperature-dependent thermodynamic equilibrium constant of one mineral phase. IAP is the product of the activity coefficients γ i and the mole fractions of the solute mineral phase x i considering their stoichiometric coefficient ν i . A positive saturation index indicates an oversaturation and a potential precipitation of the mineral. Though, if the ion activity product is smaller than the equilibrium constant, the saturation index will be negative. In this case, the solution is undersaturated with the potential to dissolve the mineral phase. Therefore, equilibrium is given at SI = 0. Debye and Hückel (1923) established an equation for non-ideal electrolyte solutions taking into account the electrostatic interaction among the ions by using the activity coefficients γ i (2). In order to fit the Debye-Hückel equation to experimental data, Robinson and Stokes (2002) extended the original equation by adding a linear concentration term ḂI.
where A and B are temperature-dependent constants, z is the charge number of the ion, I is the ionic strength, and å is the hydrated ion size. The numerator quantifies the long-range Coulomb forces acting on the ion, whereas the denominator defines the short-range interactions between the ions itself and with the solvent. The extended Debye-Hückel equation (Eq. 2) expands the former validity limit in terms of the ionic strength to I = 1.0 mol/L for mixed electrolytes. Furthermore, there are application limits given by specific temperatures and pressures. The latter is negligible at least up to a temperature of 300 °C (Helgeson 1969).
In this study, chemical speciation and saturation indices are computed with IPhreeqc 3.4.0-12927 (Parkhurst and Appelo 2013). Ion activity coefficients are based on the extended Debye-Hückel equation. The required constants A, B, and Ḃ are obtained from the commonly used LLNL (Lawrence Livermore National Laboratory) database. Saturation indices are computed for all specified minerals for a given solution. Reed and Spycher (1984) plotted these saturation indices versus temperature to investigate the equilibrium temperature of the geothermal fluid and the reservoir mineral assemblage. For validation of the multicomponent geothermometer, the method is applied on wellstudied geothermal sites in Iceland and further developed to obtain an easy-to-handle and convenient high-precision exploration tool.
Krafla is a high-temperature geothermal field located in the NE of Iceland. The geothermal system is situated in the neo-volcanic zone (Ármannsson et al. 1987). The in situ temperature measurements (Table 1) and the geochemical data of the wells (Appendix 1) were published by Guðmundsson and Arnórsson (2002). The upper 1000 m of the stratigraphy are built up by alternating layers of basaltic lavas and hyaloclastite. The latter is subglacial erupted basaltic lava, which forms hydrated breccia once it is in contact with water. Below 500 m, the hyaloclastite layers form subhorizontal reservoirs. The following 1000 m are covering basaltic intrusives, where geothermal fluids of up to 310 °C are evident (Guðmundsson and Arnórsson 2002). A more detailed stratigraphy of the field is given by Ármannsson et al. (1987). The mineralogical content of Icelandic geothermal systems and the geochemistry of the fluids are described by Arnórsson et al. (1983). Sampling methods and the geochemical analysis are given in Arnórsson et al. (2006).

Results of the analysis
The application of this multicomponent geothermometer approach requires the evaluation of the equilibrium of each solute mineral phase and, thus, the calculation of saturation indices of the considered minerals versus temperature. The saturation indices are calculated from 20 to 300 °C. The calculations of the saturation indices are processed via MATLAB. Therefore, the MulT_predict tool was developed which determines the intersection of the saturation index function for each mineral phase with the equilibrium line ( Fig. 1). Thus, the tool calculates all mineral-specific saturation indices functions throughout the temperature range by interacting with IPhreeqc. Only minerals having exactly one intersection with the equilibrium line are taken into account for the temperature determination procedure. The resulting intersection temperatures are combined in a box plot. This plot represents a first estimate of the reservoir temperature. Secondary effects perturb the chemical equilibrium of a fluid sample while migrating to the earth surface. The chemistry may change due to boiling, degassing, precipitation of phases, dilution, or mixture with shallow and low-mineralized waters as well as reequilibration with the surrounding rocks (Cooper et al. 2013;Fournier 1977;Fournier and Truesdell 1974;Pang and Reed 1998;Reed and Spycher 1984).
To determine the most vulnerable sensitive parameters for later sequential sensitivity analysis, a series of variations on system parameters (pH, redox, and steam loss) have been computed. Similarly, the concentrations of aluminium, magnesium, and iron, being major components in the minerals but only trace elements in the fluid, have been examined. The equilibrium temperature distributions for K-28 were exemplarily plotted against these parameters (Fig. 2). It is shown that the most important and vulnerable system parameter is the pH value, which is in good agreement with what is also assumed by Nitschke et al. (2017) as well as Reed and Spycher (1984). Changes have a significant impact on the solubility of mineral phases. The pH value is prone to phase segregation effects like degassing, boiling, and steam loss. Also, the pH is a temperature-dependent function, which decreases when temperature rises. In addition, regarding Fig. 2, steam loss itself is another vulnerable system parameter. Thus, possible phase segregation due to boiling has to be taken into account. The loss of steam fraction corresponds to a loss of solvent and results in the concentration of the ascending fluid. The effect of steam loss needs to be compensated by adding back the lost water. Equally, the vulnerability of trace elements is shown in Fig. 2. These constituents are particularly prone to interferences from secondary processes and measurement errors. Simultaneously, they have a high impact on the solubility product and hence, on the saturation index of the majority of reservoir minerals. Clearly, aluminium is the most vulnerable trace element. Its concentration is a crucial parameter when computing the saturation state of aluminosilicates, which represent the major mineral phases in most geothermal reservoirs (e.g. basalts, granitoids, sandstone, greywackes, etc.). Such aluminosilicate mineral assemblages contain phases like feldspars, zeolites, micas, and clay minerals. Due to the tendency of complex formation and precipitation processes (Brown 2013), the determination of accurate aluminium concentrations is prone to large errors. Furthermore, the variations of the redox potential as well as magnesium and iron concentration were tested. As it is revealed that these parameters have only marginal effects on the temperature estimations, they are not further discussed in this study. In view of the above, the in situ values of the most vulnerable sensitive parameters, pH, aluminium concentration, and steam loss, have to be reconstructed. For this optimization, a sequential sensitivity analysis for each parameter is used. This sensitivity analysis is executed by the tool. The statistically backed minimization of the temperature spread enables the back-calculation on the in situ geochemical equilibrium between the geothermal fluid and the reservoir mineral assemblage, which is the basic assumption of the method. The tool varies these parameters around the initial measured value such that a minimal temperature spread is found. In an ideal case, the equilibrium temperatures of each mineral phase of the reservoir assemblage converge to one discrete reservoir temperature. Also unknown parameters can be estimated in this manner. Thus, a geochemical foreknowledge of the geothermal system is needed to make an educated guess for the unknown parameter, which then can be estimated towards best-fit conditions. To statistically evaluate the resulting box plots, the mineral set has to remain unchanged throughout all variations. Minerals that do not equilibrate due to over-or undersaturation or have multiple intersections with the equilibrium line in any step of the sensitive analyses were discarded from further statistical processing. Changes within the set of consistent mineral phases during the sensitivity analysis would lead to false conclusions because the temperature estimations would then result from different basic conditions (i.e. different mineral sets). Concerning this, the spread of the overlaying boxes and the median differences of each neighbouring plot are matched to identify the most likely value for the sensitive parameter, indicated by the least equilibrium temperature spread. This procedure is done sequentially for the pH value, the aluminium concentration, and the percentage of steam loss. Afterwards, the best-fit values for all parameters are combined in a final temperature estimation. Nevertheless, the aim of this study is the reconstruction of reservoir temperatures, instead of the encompassing reconstruction of geochemical reservoir conditions. As a generic example, the calculation and optimization will be shown in detail for well K-28 to give an understanding of the procedure. Therefore, the geochemical data (Appendix 1) of the sample is used. The result of this first calculation is shown in Fig. 6a. For the temperature estimation from the non-specific mineral set (Appendix 2) without further optimization, a large temperature spread of 260 °C is obtained. Hence, in this study, a basalt-specific mineral set was devised to enhance the accuracy of this method. The basaltic minerals have been selected according to the mineralogical study of the Krafla reservoir rocks (Arnórsson et al. 1983). This set is extended for secondary mineral phases, occurring in geothermal reservoirs due to hydrothermal alteration processes. It is based on the stability of mineral phases at certain temperature and pressure levels which were described by Giggenbach (1981). The resulting basalt-specific mineral set (Table 2) is used to evaluate the in situ temperatures of the reservoir. After application of the multicomponent geothermometer based on the selected mineral set, the reservoir temperature estimation could be improved (Fig. 6b), yet the temperature spread still exceeds 100 °C.
As a second step, the sensitivity analysis is conducted. Firstly, the pH value is optimized. The initial pH of 9.75 is varied in increments of 0.1 towards higher acidity and basicity. The result is shown in Fig. 3 where the minimal spread of the box is reached at pH 7.85.

Table 2 Mineral phases contained in the basalt-specific mineral set devised and used in this study
Separately, the aluminium concentration is evaluated. The average aluminium concentration in the Krafla geothermal fluids is about 1.2 ppm (0.04 mmol/kg) (Guðmundsson and Arnórsson 2002). Therefore, the aluminium concentration is varied in increments of 0.006 mmol/kg. The initial aluminium concentration of the geothermal fluid composition for K-28 is 0.039 mmol/kg. In Fig. 4, an optimal concentration is also reached at 0.039 mmol/kg. Lastly, the sensitivity of fluid composition to the magnitude of steam loss is considered. The amount of steam loss is unknown, but has to be back-calculated. Therefore, pure water is virtually added back in increments of 1%. In Fig. 5, the optimum in steam loss is reached at 14%.
The final temperature estimation is then computed by combining the best-fit values for pH, aluminium concentration, and steam loss. Figure 6c displays the reduced spread of the calculated temperature after each optimization step. Simultaneously to the reduced uncertainty, the median of the temperature estimate has ascended. For validation, the concluding reservoir temperature estimations are compared to downhole temperature measurements published by Guðmundsson and Arnórsson (2002) (Table 1). Figure 7 displays the temperature box plots for the wells K-11, K-24, and for K-28 of two consecutive years. The range of the measured in situ temperatures (Table 1) is figured as an orange box. In each case, the estimated temperatures fit very well the measured borehole temperatures after the optimization procedure. Note that even very small temperature ranges are matched by the estimated temperatures (e.g. K-24 and K-28).

Discussion
The comparison of the optimized temperature estimation and the measured downhole temperatures confirms the functionality of the application. Only for K 24, the median temperature shows a minor overestimation of 1 K above the highest measured inflow temperature, though, the estimations are also located in the measured temperature range. The overall spread of each final plot after the sensitivity analyses does not exceed 7% (K-24), but is on average 3.7% of the absolute median temperature. The uncertainty throughout the validation is at maximum 2.6% of the original absolute reservoir Fig. 7 Results of wells K-11, K-24, and K-28 for three stages of the analysis. The first column displays a temperature estimation calculated based on an unspecific mineral set. The box plot in the second column represents the specified basaltic mineral set. The last box plot shows the optimized temperature estimation via the pH, aluminium concentration, and steam loss sensitivity analysis. These box plots can be compared with the range of the measured temperatures in the boreholes, given by the orange box temperature. Thus, the validated developed tool shows a significant improvement compared to uncertainties of conventional solute approaches, following in the discussion.
The calculation of the saturation indices relies on the LLNL database which is constrained to temperatures of 300 °C. Therefore, most of the geochemical modelling tools are also constrained for that p-T range. Icelandic geothermal systems have the potential to exceed these temperatures. To evaluate the validity limits, the investigations were extended to the Reykjanes geothermal system, which is located in the SW of Iceland. Furthermore, as the system is recharged by seawater, also the effects of high salinities can be assessed. Reykjanes is also situated in the neo-volcanic zone and the stratigraphy equals the Krafla structure with alternating layers of basaltic lavas and hyaloclastite in the upper part, followed by basaltic intrusives at greater depth. However, since the geothermal system is located on a peninsula, seawater infiltrates the productive horizons of the reservoir. At Krafla, the dissolved solids content is generally up to 1500 ppm. Showing sea water concentrations, the salinities at Reykjanes are up to 58 times higher (RN-23: 87.160 ppm). Óskarsson et al. (2015) published the geochemical data (Appendix 3) of the fluids from two production wells and the associated temperature logs (Table 3). In the following, the applicability of the numerical scheme will be tested at high-enthalpy geothermal fields with temperatures above 300 °C and elevated salinities.
Thus, the Debye-Hückel coefficients A, B, and Ḃ of Eq.
(2) have to be extrapolated towards higher temperatures and implemented into the database. The coefficients A and B were extrapolated by a quadratic fit, whereas Ḃ was extrapolated by a cubic fit. This extrapolation is similar to the scheme proposed by Helgeson (1969). Figure 8 shows the results of polynomial parameter estimation. The obtained values for the coefficients A and B herein are very close to the results computed by Helgeson et al. (1981). Estimations towards the critical temperature of water have to be used with care and, therefore, only exceed up to 350 °C.
These modifications of the Debye-Hückel parameters allow technically for the computation of saturation indices over an extended temperature range. To gain an overview, the saturation indices of the well RN-12 at Reykjanes were plotted to 350 °C (Fig. 9). The extrapolated coefficients follow the trend of the saturation curves. The results for the reservoir temperature estimation and the comparison against measured values are presented in Fig. 10. Despite the high sodium chloride concentrations, the tool operates thoroughly. Analogous to the methodology presented for the Krafla site, the spread of the temperature box plots is minimized and eventually matches the measured temperatures (Table 3). The spread is about 4.7% of the measured absolute reservoir temperature, while the overall temperature accuracy is at 0.5%.    (Bjarnason 2010), requiring additional gas analysis data. Afterwards, the solute geothermometers are applied. The table in Appendix 4 comprises quartz geothermometers according to Fournier and Potter (1982), Arnórsson et al. (1983), and Verma (2000); Na/K geothermometers according to Truesdell (1976), Fournier (1979, Giggenbach (1988), Arnórsson (2000), and Can (2002), as well as Na/K/Ca geothermometers according to Fournier and Truesdell (1973), Nieva and Nieva (1987), and Kharaka and Mariner (1989), and K 2 /Mg geothermometer according to Giggenbach (1988). All geothermometers were checked for their applicability in these settings. The results of the conventional geothermometers are visualized in Appendix 5 for Krafla and Appendix 6 for Reykjanes together with the results of MulT_predict. In all cases, our application targeted the measured temperature more precisely with a lower overall spread of the temperature estimation without requiring additional gas analysis data. Compared to MulT_predict (Krafla: 3.7%; Reykjanes: 4.7%), the overall temperature spread of conventional solute geothermometers is 10.5% for Krafla and 12.3% for Reykjanes.

Conclusion and outlook
This application of multicomponent geothermometry is a promising tool for reservoir temperature estimations. This specific approach comprises a devised basalt-specific mineral set and a subsequent sensitivity analysis based on a standard chemical analysis of the fluid composition without the need for a sophisticated gas analysis. The statistically robust temperature estimations of the reservoir are incorporated in a valuable tool Fig. 10 Results for wells RN-12 and RN-23 for the three stages of the analysis. The first box plot shows a temperature estimation calculated based on an unspecific mineral set. The second column displays the developed basaltic mineral set. The third box plot shows the minimized temperature estimation after the pH, aluminium concentration, and steam loss sensitivity analysis of the specified basaltic mineral set. These box plots can be compared with the range of the measured temperature in the borehole, given by the orange box for precise reservoir temperature estimation. Thus, the methodology enhances the usability, the applicability during geothermal exploration as an economically efficient tool for reservoir temperature determination.
The emphasis of this study was the validation of this optimized approach of multicomponent geothermometry. A basalt-specific mineral assemblage was devised to reduce temperature estimation uncertainties. These estimations were further improved by using a subsequent sensitivity analysis via the herein proposed MulT_predict tool. The optimization of pH, aluminium concentration, and steam loss reduces the uncertainty of the temperature estimations significantly. Hence, the back-calculations enable the reconstruction of the in situ equilibrium temperature conditions between the geothermal fluid and the reservoir mineral assemblage. This equilibrium state corresponds with the underlying geochemical assumptions of geothermometry. The approach presented here would even allow for constraining unknown input parameters. An educated guess of the parameter can be varied to reach the best-fit value. For validation, the calculated reservoir temperatures are compared against measured in situ reservoir temperatures and classic solute geothermometry. The maximum uncertainty of the temperature estimations is only 2.6% with respect to the in situ reservoir temperature. The accuracy of the results shows the efficiency and credibility of the method. This multicomponent approach benefits from its statistical robustness due to the conjunction of the saturation indices of multiple mineral phases for temperature estimations. Therefore, it can be applied to diverse geothermal sites with different fluid origins. Furthermore, high-temperature systems can be investigated by extrapolation and modification of the Debye-Hückel parameters in the LLNL database, yet resulting in reservoir temperature estimations with small variances. The method is easy to apply because of the simplicity in terms of input data. A standard water analysis is sufficient for obtaining very accurate results, which facilitates the usability especially at less explored sites where good-quality data is often missing.
Overall, the validation of the procedure was successful and improved temperature estimations via multicomponent geothermometry. In future, a broader application over different mineral sets is envisaged to expand the usability of the methodology towards other geological settings. Ystroem et al. Geotherm Energy (2020)