Understanding the Fate of H2S Injected in Basalts by Means of Time‐Domain Induced Polarization Geophysical Logging

To help meet emission standards, hydrogen sulfide (H2S) from geothermal production may be injected back into the subsurface, where basalt offers, in theory, the capacity to mineralize H2S into pyrite. Ensuring the viability of this pollution mitigation technology requires information on how much H2S is mineralized, at what rate and where. To date, monitoring efforts of field‐scale H2S reinjection have mostly occurred via mass balance calculations, typically capturing less than 5% of the injected fluid. While these studies, along with laboratory experiments and geochemical models, conclude effective H2S mineralization, their extrapolation to quantify mineralization and its persistence over time leads to considerable uncertainty. Here, a geophysical methodology, using time‐domain induced polarization (TDIP) logging in two of the injection wells (NN3 and NN4), is developed as a complementary tool to follow the fate of H2S re‐injected at Nesjavellir geothermal site (Iceland). Results show a strong chargeability increase at +40 days, interpreted as precipitation of up to 2 vol.% based on laboratory relationships. A uniform increase is observed along NN4, whereas it is localized below 450 m in NN3. Changes are more pronounced with larger electrode spacing, indicating that pyrite precipitation takes place away from the wells. Furthermore, a chargeability decrease is observed at later monitoring rounds in both wells, suggesting that pyrite is either passivated or re‐dissolved after precipitating. These results highlight that a sequence of overlapping reactive processes (pyrite precipitation, passivation, pore clogging and possibly pyrite re‐dissolution) results from H2S injection and that TDIP monitoring is sensitive to this sequence.


Introduction
Geothermal energy production emits an estimated 0.2 Mt/yr of hydrogen sulfide (H 2 S) globally, with Iceland alone emitting 30 kt/yr (Marieni et al., 2018).These anthropogenic H 2 S emissions are sourced from the geothermal production steams, which contain significant concentrations of H 2 S (Prikryl et al., 2018;Stefánsson et al., 2011).H 2 S emissions are toxic to humans and can be fatal in concentrations as low as 320 ppm (World Health Organization -Regional Office for Europe, 2000).H 2 S emissions also pose a threat to the environment, oxidizing when exposed to atmospheric oxygen to form acid rain (Greaver et al., 2012).Recent air quality regulations, both internationally and within Iceland, have limited the amount of H 2 S atmospheric emissions to 50-Iceland Ministry of the Environment, Energy and Climate, 2010).
To reduce emissions from geothermal energy production in Iceland, H 2 S can be dissolved into the geothermal wastewater and injected into the basalt subsurface (Mamrosh et al., 2014).To help meet emission standards, H 2 S from the geothermal production steam may be captured at the power plant, dissolved into geothermal wastewater, and injected into the basalt subsurface (Mamrosh et al., 2014).Geochemical and reactive transport models suggest that such a H 2 S-injection results in effective sulfide mineralization, with pyrite being the dominant sulfur bearing mineral (Bacon et al., 2014;Marieni et al., 2018;Prikryl et al., 2018;Stefánsson et al., 2011), due to the high reactivity of basalt and its divalent cations content (up to 25 wt.%Ca, Mg, Fe) (Gysi & Stefánsson, 2008;Matter et al., 2016).However, these findings are based on simulations and laboratory experiments with limited evidence from field tests (Robin et al., 2020).In practice, the H 2 S mineral storage process requires careful monitoring to identify any adverse effects of the injection, such as the acidification of the shallow groundwater system and the mobilization of toxic metals from the basaltic rocks (Cuoco et al., 2013;Delmelle et al., 2015;Flaathen & Gislason, 2007;Floor et al., 2011;Galeczka et al., 2016).
Current monitoring practices for H 2 S mineralization into pyrite are similar to those for CO 2 mineralization into carbonates (Clark et al., 2020;Snaebjörnsdóttir et al., 2017).They consist mainly of quantification by (a) mass balance calculations using tracer tests (Matter et al., 2016), (b) transport models (Ratouis et al., 2022), and (c) steady-state "reaction path" geochemical models based on chemical monitoring.Chemical and mineral analyses on precipitates formed on metallic surfaces within monitoring wells also bring qualitative mineralization evidence (Matter et al., 2016).More recently, the triple sulfur isotope systematics (32, 33, and 34S) of geothermal fluids were investigated as a possible tool to trace H 2 S sequestration by sulfide mineralization in the geothermal reservoir (Robin et al., 2020).One major pitfall of these monitoring practices is that they typically capture less than 5% of the flow and extrapolate the results (Matter et al., 2016), whereas the remaining 95% may include undetected upward migration through fractures and springs.In addition, a key assumption behind the quantification is that the amount of dissolved gas not captured by the measurement downstream, compared to expectations, has mineralized and that mineralization is permanent (Gunnarsson et al., 2018;Matter et al., 2016).While physical verification of these assumptions, for example, through post-injection core-drilling, has not yet been possible (Carbfix, 2022), other studies focusing on CO 2 indicate that additional processes to carbonate mineralization can happen, namely the formation of unstable organic carbon due to microbial activity (Trias et al., 2017).Therefore, quantifying mineralization with these methods, whether sulfide or carbonate, leads to considerable uncertainty (White et al., 2020).
Geo-electrical methods are particularly suitable for pyrite detection.Underground massive metallic ore deposits, as well as veinlets and disseminated sulfides, have been detected and discriminated with induced polarization (IP) field measurements (Börner et al., 2018;Pelton et al., 1978;Placencia-Gómez, 2015).Petrophysical laboratory measurements also indicate that the presence of pyrite can be detected in altered volcanic rocks and quantified to some extent, with induced polarization (Lévy, Gibert, et al., 2019).The polarizability, also called chargeability, describes the amplitude of electrical polarization in rocks, which is particularly strong at interfaces between fluid and semi-conductors, such as pyrite or magnetite (Abdulsamad et al., 2017;Bücker et al., 2018;Gurin et al., 2015).The reorganization of charges at the interface between pore fluid and metallic particles creates local electrical "displacement" currents, which add up until a maximum voltage is reached.The time-delay between maximum voltage and maximum current, also known as the phase-angle for sinusoidal signals, is caused by the polarization (Bücker et al., 2018, ;Chelidze & Gueguen, 1999;Olhoeft, 1985).On the other hand, the conductivity quantifies the amplitude of the electric conduction phenomenon in rocks, which is in particular enhanced by the presence of smectite, as well as high porosity, salinity, and temperature (Flóvenz et al., 2005;Kristinsdóttir et al., 2010;Lévy et al., 2018;Waxman & Smits, 1968).Geo-electrical field investigations carried out at a geothermal site in North-East Iceland confirmed that the conductivity and IP signals of naturally present pyrite, smectite, and magnetite can be picked up down to 200 m depth (Lévy, Maurya, et al., 2019).
While many applications of geophysics focus on structural characterization using static geophysical data (Binley et al., 2015), time-lapse geophysics provides a promising avenue to characterize the spatial distribution and temporal evolution of transport and reaction processes (Hermans et al., 2023).So far, most of the field applications of time-lapse geophysics have focused on characterizing flow and transport processes: groundwatersurface water interaction (Slater et al., 2010), water absorption by roots (Jayawickreme et al., 2010;Mary et al., 2020), migration of methane plume (Steelman et al., 2017), migration of injected remediation agent at contaminated site (Flores Orozco et al., 2015;Lévy et al., 2024).Monitoring reactive processes with time-lapse geophysics remains challenging and a very limited number of field cases exist (see e.g., reviews by Binley et al. (2015) and Hermans et al. (2023)).An important field demonstration was carried out by Williams et al. (2009), where time-lapse IP was used to infer subsurface accumulation of iron-sulfide minerals due to microbially-mediated iron and sulfate reduction, down to 10 m depths.While this study illustrates the sensitivity of time-lapse IP surveys to iron-sulfide precipitation and evolution over time, the generation of IP anomalies accompanying the injection of H 2 S in basaltic geothermal reservoirs is poorly understood.In addition, the geophysical survey type investigated by Williams et al. (2009) is irrelevant to investigate the depths of interest in the context of H 2 S injection at geothermal sites.
In January 2021, the geothermal power company ON Power started a small-scale injection of H 2 S and CO 2 dissolved in seal water from liquid-ring vacuum pumps at Nesjavellir, South-West Iceland.This small-scale injection runs alongside larger-scale injection of separated geothermal fluid from the power plant that has been part of standard operations since 1990, shortly after the power plant was installed.The injection takes place in eight shallow wells ranging from 311 to 660 m deep.We use this injection experiment to investigate (a) whether time-lapse geophysics is able to detect changes related to H 2 S injection and how significant these changes are compared to noise and (b) what novel insights the time-lapse IP method can offer regarding potential complexities in the sustainability of pyrite mineralization not previously considered.We also evaluate (c) the applicability of the existing Debye decomposition approach to calculate relevant parameters for converting the changes of IP signal into volumes of pyrite formed, in a context where power line noise affects data quality, in particular by distorting the shape of IP decays.As surface geophysical measurements lose resolution with depth, wireline logging was considered a relevant intermediate scale between laboratory and surface geophysical measurements to assess pyrite precipitation in the vicinity of the injection wells.The focus of our study is thus to evaluate the potential of time-lapse IP through state-of-the-art wireline logging, as a complement to geochemical monitoring (Clark et al., 2020;Snaebjörnsdóttir et al., 2017), to monitor H 2 S mineralization and its evolution over time, down to 550 m depth.

Electrical Response of Pyrite With Frequency-Domain IP (FDIP)
In an electrical context, geological formations can be thought of as a combination of resistors and capacitors (Lévy, 2019).Ions in pore water and those connected to clay minerals are charge carriers in aquifers.Conductive water and clay minerals can be represented by conductors with high conductivity, as opposed to freshwater and pure quartz sand or gravel, which have lower conductivity.Another type of good electrical conductor is common in geological media: semi-conducting minerals (e.g., pyrite or magnetite).
Conduction in most semi-conducting minerals is made possible by carrier traps, due mainly to (a) vacancies (missing atoms), (b) substitutional impurities (replacement by an atom of another element, such as a trace element).This allows generating an electron-hole pair, by changing an electron from an orbital with energy near the top of the valence band to one with energy near the bottom of the conduction band (Revil, Florsch, & Mao, 2015;Shuey, 1975).Indeed, ions of transition elements, such as Fe or Cu, have partially filled 3d electron shells, meaning that some of the 3d orbitals contain electrons only part of the time (Shuey, 1975).In this context, electrons in the normally empty conduction band or holes in the normally full valence band are free charge carriers (Revil, Florsch, & Mao, 2015;Shuey, 1975).In the case of pyrite, the overlap of atomic orbitals across atoms, associated to covalent bonds in the crystal, allows a net motion of electrons and holes within the crystal, from atom to atom (Shuey, 1975).Under an alternating electrical field, charges move back and forth.For the free charge carriers in semi-conducting minerals, the limit on the amplitude of charge displacement is on the order of the mineral size, that is very large compared to the displacement of bound charges, for example, in quartz (Shuey, 1975).The motion of charge carriers throughout the crystal enables electro-diffusion of electrons and holes over the whole crystal length.This results in strong polarization effects on the margins of the crystals, and thus a strong chargeability at low frequency (Abdulsamad et al., 2017;Revil et al., 2018;Shuey, 1975;Wong, 1979).In geological media semi-conducting minerals are often disseminated and thus mainly responsible for a strong chargeability, that is, induced polarization, also measured by the phase shift between the current and the voltage in frequency-domain.Beyond approximately 20%-30% in volume, semi-conducting minerals may be interconnected over sufficiently long distances (e.g., in veinlets) to also cause a strong conduction (Shuey, 1975;Slater et al., 2005;Slater et al., 2006;Revil, Abdel Aal, et al., 2015).According to Gurin et al. (2019), mechanisms related to the chemical interface between semi-conducting minerals and pore electrolyte also contribute to strong induced polarization.Overall, the question on which phenomena dominate in natural geological media is still explored (Abdulsamad et al., 2017;Izumoto, 2023;Shuey, 1975).
Figure 1 illustrates the large polarization, represented by the phase angle parameter, of volcanic samples containing pyrite, as opposed to a volcanic sample of similar geology and full of smectite but without pyrite; the latter being more conductive but much less polarizable, as described in details by Lévy, Gibert, et al. (2019) and Lévy, Weller, and Gibert (2019).It is also visible in Figure 1 that the grain size distribution of metallic particles affects the frequency at which the phase angle peak occurs.The laboratory results presented in these two studies are obtained on Icelandic volcanic samples with various types of semi-conducting minerals (pyrite, pyrrhotite, magnetite), grain size and quantity, as well as variable porosity and clay type and content.Since volcanic rocks at Nesjavellir share characteristics and features with commonly encountered Icelandic volcanic rocks (e.g., high amount of volcanic glass, basalt chemistry, type of porosity), laboratory results obtained on samples at another Icelandic site are used as a basis for the present field study.
Several models have been developed to describe how ions reversibly accumulate at the interface between pore fluid and metallic particles, causing polarization.One of the most widespread phenomenological relaxation models to interpret FDIP data is the Pelton model (Pelton et al., 1978;Tarasov & Titov, 2013;Weller & Slater, 2022), given by Equation 1. (1) with four fitting parameters: the DC resistivity ρ DC , the chargeability m, the Pelton time constant τ and the socalled Cole-Cole exponent c.
However, the Pelton model, presented in Equation 1, is only valid for spectra measured on samples with a single grain size fraction.Variations of this model for two grain size fractions exist but the fitting parameters are not directly comparable to each other, and even the two grain size fractions model may not represent the complexity of naturally occurring pyrite crystals in volcanic rocks (Lévy, Gibert, et al., 2019).
A more general approach to analyzing FDIP data is the Debye decomposition (Nordsiek & Weller, 2008).With this approach, the frequency-dependent complex electrical resistivity ρ*( f ) describes an IP spectrum that can be fitted by a superposition of N Debye terms (Pelton models with c = 1), as shown in Equation 2.
where m k and τ k are pairs of partial chargeability and relaxation time of a single Debye model.The function m k = f (τ k ) is also referred to as relaxation time distribution (RTD).The N chargeability values m k , which are related to a set of pre-defined relaxation times τ k , and the value of DC resistivity ρ DC are the result of a fitting procedure described by Nordsiek and Weller (2008).Two parameters can be integrated from the Debye decomposition to simplify the interpretations: the total chargeability m tot (Equation 3) and the mean relaxation time and τ mean (Equation 4).
A relationship linking the chargeability to pyrite volume in the absence of background polarization was developed by Revil, Florsch, and Mao (2015), based on effective medium theory(Equation 5).
where p v is the volume fraction of pyrite (or other metallic particles) and m is the chargeability in V/V (dimensionless).In their study, the chargeability parameter is fitted to frequency-domain IP (FDIP) data by a Pelton model (Cole & Cole, 1941;Pelton et al., 1978).According to Martin and Weller (2023), the chargeability m fitted by a Pelton model (Pelton et al., 1978) and m tot fitted by Debye decomposition (Nordsiek & Weller, 2008) are similar for sand-pyrite mixtures.
While the influence of temperature on the polarization effect is still poorly understood, recent laboratory experiments indicate that the phase angle in the presence of metallic particles is temperature independent in the range 5-50°C (Revil et al., 2018).

Frequency-Domain and Time-Domain Induced Polarization (FDIP and TDIP)
In the field, the Time Domain Induced Polarization (TDIP) method is often preferred to FDIP, primarily due to its easier and faster implementation.TDIP consists of injecting a direct current pulse of finite duration (DC) through two current electrodes (A and B) and recording the resulting voltage at two potential electrodes (M and N) with a fine-enough sampling rate to capture the exponential charge and discharge of the voltage, during the on-time and off-time, respectively (Sumner, 1976;Telford et al., 1990).Further details on the TDIP methods can be found in Binley and Slater (2020).
In general, voltage discharge measured by TDIP is expected to have a stretched shape, which describes the relaxation behavior in a similar manner as Cole-Cole and other related FDIP models (Alvarez et al., 1991;Cole & Cole, 1941;Davidson & Cole, 1951;Havriliak & Negami, 1966).More details are described in Appendix B.
While quantitatively connecting the TDIP and FDIP responses of rocks containing metallic particles is not trivial (Alvarez et al., 1991), the Debye decomposition is a common method for analyzing both TDIP and FDIP data (Nordsiek & Weller, 2008;Tarasov & Titov, 2007).An equivalence was found by Martin et al. (2021) between the Debye decomposition of FDIP and TDIP data, provided that the current injection time in TDIP was long enough (ideally 64 s, but decent comparisons were found beyond 2 s).With 2-s current injection, as is the case in this study, the coefficients from the Debye decomposition of TDIP data are slightly underestimated compared to FDIP data.Despite this limit, this approach was considered the most relevant one.As opposed to the integral chargeability, which tends to smoothen out noise in the data, the Debye decomposition tends to amplify the noise and, thus, requires careful data processing (Martin et al., 2021).

H 2 S Production and Re-Injection at Nesjavellir Geothermal Field, Iceland
The methodology developed in this study for monitoring H 2 S injection was tested at the Nesjavellir geothermal field (SW Iceland), where a geothermal power station has been operated since 1990 (Figure 2).Nesjavellir is a high-temperature geothermal field (>200°C above 1 km depth) characterized by a relatively high permeability (Zakharova & Spichak, 2012).Nesjavellir is located just north of the Hengill central volcano, which is situated on a ridge-ridge-transform triple junction (Foulger, 1988).The area is comprised of mainly late Quaternary hyaloclastites and post-glacial age hyaloclastites and basalt flows (Árnason et al., 1969;Foulger & Toomey, 1989).This study focuses on shallow depths northeast of the geothermal production area and above the low-permeability clay cap layer (approximately 500 m depth) (Franzson & Gunnlaugsson, 2020;Gómez-Díaz et al., 2022;Gunnarsdottir et al., 2020;Schiffman & Fridleifsson, 1991).Intrusions are rare at these shallow depths, and none have been identified in the upper 300 m (Franzson, 1988).
The Nesjavellir power station currently produces 120 MWe of electricity and 300 MWth of thermal energy for district heating.The hydrothermal production fluid (260-300°C) is sourced at 1,000-1,500 m depth (Snaebjörnsdóttir et al., 2020).Following a controlled pressure-decrease, the hydrothermal fluid is divided into (a) steam and (b) geothermal separated water at 192°C.Electricity generation at Nesjavellir uses 240 kg/s of steam to power four turbines.The steam is then condensed, using 2,000 l/s of cold groundwater (5-7°C), to form condensate wastewater, composed primarily of distilled water.In parallel, heat exchangers use the 192°C geothermal separated water to heat cold groundwater to 87°C, which is then pumped to the capital Reykjavík, 27 km away.The separated water contains the dissolved solids from the hydrothermal fluid and thus has a high potential for chemical pollution and corrosion.Reaction path modeling from the literature predicts that basaltic glass dissolution is the main source of iron over short time durations, according for example, to Equation 6, which corresponds to stoichiometrically-defined basalt dissolution in the Carbfix database (Heřmanská et al., 2022;Stefánsson et al., 2011).Magnetite dissolution contributes less iron due to the slower dissolution kinetics and more limited amount.
Pyrite precipitation is expected to be rapid (Appelo & Postma, 2005;Prikryl et al., 2018;Stefánsson et al., 2011), following a 2-step reaction, presented in Equation 7 (Appelo & Postma, 2005). 2FeOOH An overview of the area, including roads, boreholes, and buried infrastructure, is presented in Figure 2. Information on the wells and the injection are given in Table 1.Temperature logs in NN3 and NN4 measured the days of TDIP logging are presented in Figures C1 in Appendix C, together with caliper logs measured after drilling (Figure C2).Further details on the wells can be found in Gómez-Díaz et al. ( 2022) and Hafstað (2003).

Field Monitoring With QL40-IP Logging Tool
The QL40-ELOG/IP (Mount Sopris QL40-IP, 2020) logging tool measures the electrical resistivity and timedomain IP (TDIP) response with electrodes made of stainless steel.Stainless steel electrodes are commonly used in TDIP field studies, as the effect of electrode polarization is considered to have negligible impact on data quality for current injection time lasting 2 s or less (Dahlin et al., 2002).The tool has a diameter of 43 mm and uses the "normal" electrode configuration with 16-and 64″ electrode spacings, defined as the distance between the current electrode (A) and the potential electrode (M) (Helander, 1983).The resistivity is obtained after multiplication of the resistance by the geometric factor 4πd AM .The current generator, built into the tool, sends current at electrode A in all directions and the shielding of the logging wireline serves as the current sink.The reference potential electrode (N) is located on top of an 8 m-long isolation bridle, such that the distance between the potential electrodes (M and N) is considered infinite in comparison to the much smaller electrode spacing AM.Moreover, grounding of the logging unit (truck) is achieved by clamping the truck ground onto the casing for security reasons (Advanced Logic Technology, 2021).The 64″ measurements are sensitive to polarization effects over a larger measurement volume, deeper into the formation, compared to the 16″.The radial investigation characteristic, as defined by Roy and Dhar (1971), considers electrostatic potentials of individual cylindrical shells of varying radii integrated to obtain the total contribution to the measured signal.According to the radial investigation characteristics formula by Roy and Dhar (1971), that utilizes the injected current and measured resistivity values to define the relationship between signal and distance, 75% of the measured signal with the 64″ spacing is caused by the cylindrical volume comprised in a 2.5-m radius from the borehole central line.As a comparison, 75% of the measured signal with the 16″ spacing comes from within 0.5 m, which is about twice the boreholes' mean radius (Table 1).
Here, measurements took place every 25 cm, with a logging speed of 1.8 m/min.
A total of eight measuring rounds took place but only four measuring rounds are presented in this study for NN-3 and NN-4 due to mainly one problem encountered on-site, related to the inability to inject current higher than 200 mA, happening randomly at certain places and dates (Table 2).
Traditional applications of the QL40-IP tool, for example, in the mining industry, are typically fulfilled with 250 or 500 msec current injections (on-time) and a fixed number of 400 sample points in the whole cycle T on,+ ; T off,+ ; T on, ; T off, .Two important differences between traditional applications and the present study led to the development of a new processing board.First, Martin et al. (2021) show that FDIP and TDIP data measured on the same samples decently overlap when current injections last at least 2 s during TDIP data acquisition.This overlap is key to quantifying pyrite volumes by Debye decomposition since these empirical relationships have only been established with FDIP data.Second, the 132 kV buried power line connected to the Nesjavellir power plant (Figure 2) is responsible for a strong 50 Hz signal in the TDIP data, as revealed by initial tests in NN-4 prior to this monitoring study (December 2019).These tests were carried out with the original QL40-IP tool, using 250-msec and 500-msec current injections, corresponding to full cycles of 1 and 2 s, respectively, and thus sampling rates of 400 and 200 Hz (400 sample points over the whole cycle).Using 2-s current injections, the full cycle lasts 8 s, which leads to a sampling rate of 50 Hz with 400 sample points and aliases the 50 Hz noise.Further investigations also showed that both the amplitude and exact frequency of the "50 Hz" noise varied over time.Therefore, it was clear that TDIP data acquisition with 2-s current injection required a larger number of sample points to model the noise and remove it before data analysis.This was achieved by increasing the memory of the processing board and led to an updated version of the QL40-IP tool used in the rest of this study.
From September 2020 (baseline) and on, the following IP acquisition settings were used: square-wave current injection (T on,+ ; T off,+ ; T on, ; T off, for 2 s each, 8 s total), 450 Hz sampling rate of the voltage during the whole cycle (i.e., 3,600 sample points), and wireline speed of 1.8 m/min.The 50 Hz noise is further illustrated in Figure 3 (first and second columns).
In parallel to the TDIP monitoring, the changes in electrical resistivity were also monitored to evaluate the possible contribution of processes other than pyrite precipitation to the changes in the IP response.Indeed, changes in borehole fluid composition, temperature, and precipitation of clay minerals are reflected by changes in resistivity.Similarly, the formation of connected clusters of pyrite particles that would lose their polarization properties while becoming more conductive would also cause changes in the resistivity logs.

Data Processing From Full Waveform Voltage Signal to Polarizability Decays
The apparent chargeability M app (t), also called polarizability and noted η(t) in Martin et al. (2021) and Tarasov and Titov (2007), corresponds to the ratio between the decaying voltage during the off-time, V decay (t), and the maximum voltage reached during the on-time, V DC .It is calculated by stacking decays related to positive and negative current polarities (Equation 8).A series of processing steps are carried out to clean the full waveform signal, optimize the signal-to-noise ratio, and maintain three time-decades of signal for the polarizability.The term polarizability is used in the rest of the paper to refer to the apparent chargeability defined in Equation 8 in order to avoid confusion with the integral chargeability or the total chargeability parameter calculated from Debye decomposition.First, a time-invariant self-potential (SP) voltage affects both V(t) and V DC in the positive and negative decays.
The SP offset at each depth and each date is calculated by averaging the positive and negative DC voltages, V DC,+ and V DC, and removed from V(t) and V DC .
Second, harmonics of 50 Hz noise, due to a 132 kV buried power line, affect the voltage full waveform measurements, as mentioned in the previous section and illustrated in Figure 3.A power line noise model, based on the algorithm proposed by Larsen et al. (2022) and applied in a similar manner by Olsson et al. (2016), is used to subtract the harmonic noise from the signal.As shown by Larsen et al. (2022), the fundamental frequency varies over time, rapidly oscillating between 49.9 and 50.1 Hz, and must first be determined.The 2-s decaying voltage is thus decomposed into 200 msec segments, where the frequency is fitted in every segment.Once the exact frequency is determined, a model of the background noise, including harmonics up to Nyquist frequency, is removed from the signal.The denoised voltage after applying this 50 Hz filter is presented in Figure 3 (panels in the second column) as the red curve on top of the FW voltage (in black) during the positive off-time section of the full cycle.
Third, the voltage decays are re-gated into log-spaced time windows.This step transforms the full waveform voltage decay signal with 900 linearly spaced data points into a gated voltage signal with 37 logarithmically spaced gates.More accurately, the re-gated voltage consists of 10 gates in the first 22 msec (2.21 msec width, i.e., the original sampling step) and then exponentially larger gate widths (from 4.42 to 141.44 msec), over which an exponentially increasing number of data points (from 2 to 64) are averaged.Given the strong exponentially decaying shape in the early times, it appeared best to keep all data points in the first 22 msec.Since the original data points are kept as the 10 first gates, it corresponds to a "square-gate" with one data point averaged.To be consistent, the square gating strategy is kept throughout the decay.This is especially relevant because the signal can be locally approximated by linear functions in the later times.The gating strategy is further illustrated in Figure A1 in Appendix A. At some depths, the first gate (centered at 1.1 msec) is affected by electromagnetic (EM) coupling (see, for example, at +263 m in Figure 3).Therefore, it was discarded from all the re-gated discharge curves, and the first effective gate is at 3.32 msec: visual inspection of all the decays suggests that EM coupling is negligible at this gate, which is consistent with observations by Olsson et al. (2016).
Fourth, the positive and negative voltage decays are divided by the DC voltage and stacked to obtain the polarizability in mV/V (Equation 8).The stacking allows removing extra noise with the same sign in the positive and negative decays.These four processing steps eventually provide the processed polarizability curves over almost three decades (in the range of 3-2,000 msec), as presented in Figure 3 with a linear scale (third column) and log-log scale (fourth column).The linear scale allows better comparison to the full waveform signal presented in the same figure, while the log-log scale is used in the rest of this study.
Fifth, polarizability decays are averaged over 2-m thick sections (Figure 4), corresponding to 8 decays averaged together (0.25 m spatial resolution).This is done to facilitate the comparison between measuring rounds since the exact depth of each single measurement varies at different dates.It also allows smoothing out short-wavelength signal oscillations that would prevent meaningful comparison between measuring rounds.The standard deviation corresponding to this average also captures noisy areas where the signal varies a lot over short distances.The averaging over 2-m sections also takes into consideration the fact that logging measurements at one depth are, in reality, sensitive to a sphere around this given depth.Finally, averaging over 2-m sections helps to identify and mitigate remaining harmonic noise.
Some polarizability curves remain noisy after the processing (non-monotoneously or non-exponentially decaying), for example, at 263 m depth and +40 days shown in Figure 3 and at 382-384m and +40 days shown in Figure 4. Noisy processed polarizability curves are further addressed in the next subsection and discussed at the end of the paper.

Estimations of Pyrite Precipitation From Polarizability Decays
A simple way to represent TDIP data (polarizability decays) as a function of depth is to integrate the decays to obtain the so-called "integral chargeability", M int (Mao et al., 2016;Telford et al., 1990).It represents the area enclosed by the discharge curve, V(t), and its zero asymptote, in a given time-window, [t1: t2], divided by the primary voltage V DC and the time-window width, as described by Equation 9 (Bertin & Loeb, 1976;Sumner, 1976).An error on M int , Err(M int ), is calculated following Equation 10 at every date d.
where std(η i ) is the standard deviation of the polarizability of the i th gate within the averaged 2-m thick section, n is the number of time gates in the polarizability decays and η is the average polarizability over the n gates.
The integral chargeability presents the advantage of smoothing out the noise.It is calculated here in the time interval ranging from the second gate (3.32 msec) to the last gate (1,926 msec).M int is used in this study as a qualitative assessment tool of the change in polarization over time at a given depth, the trend corresponding to pyrite formation or dissolution.The error bar of the relative difference Δ rel is calculated with Equation 11for the integral chargeability M int and the resistivity ρ.
where d and d 0 correspond to the current monitoring date and the baseline date, respectively, Err(ρ(d)) is the standard deviation within the averaged 2-m thick section, and Err(M int (d)) is calculated as per Equation 10.
The integral chargeability calculated from TDIP data cannot be converted to pyrite volume as M int is not equivalent to the "Cole-Cole" or "Pelton" m in Equation 5. Several approaches have been developed to relate TDIP and FDIP data, considering there is no strict equivalence between Pelton models in frequency-domain and stretched exponential functions in time-domain (Alvarez et al., 1991).As described in Section 2, there is an equivalence between the Debye decomposition calculated with FDIP and TDIP data, with a slight underestimation of m tot from TDIP data when the current on-time is 2 s (Martin et al., 2021).Therefore, the conversion to pyrite volume fraction is carried out here by using Debye decomposition.The polarizability M app = η(t) is fitted by a non-negative linear least-square procedure, where the fitting function is given in Equation 12and the optimization problem is expressed in Equation 13.
where x is a vector with the m k model parameters, d is a vector with the measured data η(t i ) and C is the linear multiplier defined as C i,k = e t i τ k , where t i are the center times of the p gates (i = 1..p) and τ k are the pre-defined N relaxation times (k = 1..N).Here we define the τ k,k = 1..N as a set of N = 300 logarithmically-spaced values of discrete relaxation times in the range 3-2000, covering approximately three decades, following the methodology proposed by Nordsiek and Weller (2008) and also applied by for example, Robinson et al. (2018) and Z. Zhang et al. (2017).The range 3-2,000 msec considered for the relaxation times corresponds to the measuring window.During the procedure, most of the coefficients m k are kept to 0 due to the non-negative constraints.The fit eventually results in 10-20 coefficients being non-zeros.
The total chargeability m tot is calculated by summing all the m k coefficients (Equation 3).The difference of m tot at different measuring rounds, Δm tot , is converted into a difference in pyrite volume fraction p v , Δp v , following Equation 14, which is based on Equation 5.
where d and d 0 are the monitoring date and the baseline, respectively.Only the difference is converted, not the absolute pyrite volume fraction, in order to subtract the contribution of the background polarization that may be partly due to iron-oxides or other polarizable material present before the start of H 2 S injection.The assumption that only pyrite is responsible for the changes in the IP signal is further discussed in Section 5.1.
As opposed to the integral chargeability, the Debye decomposition seems to act as a noise amplifier, similar to observations by Martin et al. (2021).As mentioned previously, some polarizability curves remained noisy despite the numerous processing steps, with a resulting non-decaying shape.Such curves should not be processed with Debye decomposition.An additional processing step is introduced to ensure that only exponentially-decaying curves are processed with Debye decomposition.The criterion for ensuring this shape is based on the Kohlrausch-Williams-Watts relaxation model for TDIP data (Alvarez et al., 1991), presented in Appendix B. In practice, the polarizability curves η(t) are fitted with a stretched exponential function, presented in Equation 15, before applying the optimization procedure described in Equations 12 and 13.
The fitting parameters α, β, and τ KWW , further described in Appendix B, are not used for interpretation in this study.Instead, the stretched exponential fit serves two purposes: (a) calculating a deviation to a stretched exponential to remove outliers polarizability curves where R 2 < 0.998 and (b) smoothing out the polarizability decays that are reasonably deviating to ensure the Debye decomposition is not fitting noise.In order to confirm that this fitting procedure is not introducing a significant bias, the m tot and τ mean resulting from Debye decomposition before and after the fitting procedure are compared in Figures B1 and B2 in Appendix B. While other types of curves may reflect IP effects, with for example, negative polarizability values having a physical explanation in the context of dipole-dipole acquisition with sharp layer boundaries (Fiandaca et al., 2022), the polarizability curves discarded here did not belong to these "heterodox transients," but rather had a distorted shape, similar to the curve in NN3 at +40 days at 382-384 m depth, presented in Figure 4.
An error is calculated on the resulting values of m tot (Equation 3), taking into account both the residuals to the stretched exponential fit and the standard deviation within the averaged 2-m sections, to help discriminate between variations that fall within an uncertainty interval from significant variations over time.First, an error is calculated for each gate of the polarizability curve as the sum of the 2-m standard deviation (corresponding to the variation within the averaged 2-m section, see also Equation 10) and the residuals to the stretched exponential fit (Equation 16): if the residual is positive, it is added to the positive error bar Err + , if the residual is negative, it is added to the negative error bar Err .
where η i,fit is the value of the stretched exponential fit of the polarizability at the ith gate, and Err + and Err are the positive and negative error bars on η i, fit at each gate.
Similarly to the error on M int , a relative error is calculated for each decay using the root mean square of the error on all the gates divided by the average polarizability η.This relative error is then applied to the total chargeability m tot , then to Δm tot (change from the baseline), and finally to the pyrite volume fraction change Δp v , following Equation 17.

Baseline Resistivity and Integral Chargeability
In order to evaluate the amplitude and meaning of the changes in the following sections, the resistivity and integral chargeability measured in the baseline are first presented in Figure 5.The average resistivity of NN3 and NN4 is about 30 Ωm, with a clear threshold at 450 m in NN3, separating the well into two sections, above and below 30 Ωm.The overall weak contrast between the fluid conductivity (10-15 Ωm, see Table 1) and the rock conductivity (30 Ωm) indicates that the desired current can be injected into the formation at most depths.Both in NN3 and NN4, a few intervals have resistivity values significantly higher than 30 Ωm, up to 200 Ωm.These high-resistivity intervals coincide with low-chargeability intervals (gray rectangles in Figure 5).The average integral chargeability is around 20-30 mV/V, with generally higher values in NN4.
Intervals of higher resistivity correspond to more compact, less porous, lithological layers, where less hydrothermal fluid flow and thus alteration is expected (Lévy et al., 2018).Therefore, the relatively high chargeability in the baseline, as well as the coincidence of locally high resistivity with locally low chargeability, suggest that hydrothermal alteration had been occurring in the wells prior to the start of H 2 S injection, especially in intervals of low resistivity.Small amounts of magnetite are expected in unaltered zones, typically 2-4 wt.% (Lévy et al., 2018) and may contribute to the chargeability in high resistivity layers (Lévy, Gibert, et al., 2019;Peshtani et al., 2022).
Alteration minerals, such as pyrite and smectite, are especially expected in permeable, fractured layers (Gudmundsson et al., 2010;Lévy et al., 2020;Lévy, Gibert, et al., 2019) and may be responsible for the higher chargeability in low resistivity intervals, particularly below 450 m in NN3 and below 375 m in NN4.Note that the relative contribution of polarization (quadrature conduction, representing displacement currents) to the overall conduction (in-phase and quadrature) reduces when the volumetric smectite content in volcanic rocks exceeds 10% (Lévy, Gibert, et al., 2019;Lévy, Weller, & Gibert, 2019;Mendieta et al., 2021;Telford et al., 1990), as illustrated in Figure 1.Based on these observations, the interval below 450 m in NN3 is where extensive pyrite precipitation is expected following H 2 S injection, which is also consistent with the observation that most of the "natural" hydrothermal alteration happens below 450 m in the area (Helgadóttir, 2021).
In both wells, higher chargeability is observed with the 64″ spacing compared to the 16″.Considering that the fluid is not chargeable in this frequency-range, this illustrates that the 64″ images more volume of the chargeable properties of the rock formation, and thus more changes are also expected with the 64″ spacing.On the other hand, similar resistivities are observed with both spacings, which reflects the fact that the fluid and the formation have resistivities in the same range.Based on fluid conductivity logging, the borehole fluid resistivity is around 10 Ωm in both wells (Table 1), that is, only three times lower than the average 30 Ωm measured for the formation.

Changes in Resistivity and Integral Chargeability Over Time
Logging results at +40 days (NN3 and NN4), +270 days (NN3), +380 days (NN4), and +540 days (NN3 and NN4) are presented as relative differences (in percentage), compared to the baseline in Figure 6.The average value within a 2-m thick interval is taken for each round to smooth out rapid vertical variations.Error bars, calculated from the standard deviation of this 2-m average (Equation 11), allow discriminating significant changes over time from changes within the local variability.The resistivity after baseline is further corrected to remove the temperature effect, following Equation 18 and using temperature logs shown in Figure C1 in order to assess resistivity changes reflect changes in porosity or clay content.
where T and T 0 are the current temperature and the reference temperature (corresponding to the baseline temperature), respectively.And α T = 0.02°C 1 is based on Arps (1953).The effect of temperature on IP processes is still not completely understood but laboratory investigations on sandstone, as well as samples containing semiconducting minerals report no or very limited influence of temperature on the chargeability (Revil et al., 2018;Zisser et al., 2010).Influence of temperature on the relaxation time may be described by similar laws as Equation 18 which corresponds to limited variations, compared to other influencing factors with high uncertainty such grain size, grain shape, grain orientation, porosity and pore size.Therefore, the influence of temperature on IP parameters is disregarded here.
It can first be observed that changes depend on the electrode spacing for both the resistivity ρ and the integral chargeability M int : relative differences are larger with the 64" spacing, indicating that more changes are happening in the rock formation than near the borehole wall (Figure 6).Resistivity variations are in many places within error bars, and thus, less focus is given here to interpret these variations.
Significant changes in the integral chargeability are observed.In NN4, M int increases at +40 days by 10%-20% with the 64″ spacing and 5%-10% with the 16″ spacing (Figure 6).At the next monitoring date, that is, +380 days,  M int decreases back to the baseline in most places with both spacings (note that data above 330 m were discarded for this date due to low current injected).At the last monitoring date, +540 days, M int remains at baseline values down to 300 m, but the 64″ signal decreases below baseline values between 300 m and the bottom of the well (380 m).The 16″ signal mostly remains at baseline values all along the well.A similar trend is observed in the bottom part of NN3, below 450 m: first M int increases at +40 days by up to 50% with the 64″ and up to 20% with the 16″ spacing.Then, M int decreases for both spacings but remains above baseline values everywhere, unlike in NN4.At +540 days, below 450 m, slightly higher M int values are observed compared to +270 days for the 64″ spacing.The 16″ signal follows the same trends everywhere, with lower amplitude.The maximum M int increase is smaller in NN4 than in NN3, but the increase is more uniformly spread along the borehole in NN4.In particular, a significant increase is observed at +40 days in the whole well, that is, between 200 and 380 m.In NN3, the most important changes are located between 450 and 550 m, with only minimal changes down to 450 m.
A consistent observation links the 16-and 64" spacing in NN3 and NN4: M int does not follow a monotonous increase as would have been expected if pyrite had been precipitating continuously or even if pyrite had precipitated up to a certain amount and remained in place.To illustrate further the non-monotonous pattern observed in the four panels showing M int variations in Figure 6, the evolution of M int with time in certain depth ranges is presented in Figure 7 for NN3 and NN4, using only the 64″ spacing (larger variations) and showing absolute differences in mV/V.
Resistivity changes in NN3 with the 64″ spacing show a spatial correlation with M int changes: a decrease of resistivity is observed in the upper part (above 450 m), and an increase is observed in the bottom part.As opposed to the time-evolution observed for M int , the resistivity increase at the bottom of NN3 seems monotonous (Figure 6).In NN4, a resistivity decrease is observed at +380 days above 300 m, and an increase is observed between 300 and 350 m.However, an almost uniform resistivity decrease is observed at +540 days, both with the 64″ and 16″ spacing, which indicates that the fluid conductivity may have changed due to a change in fluid composition.Since the fluid conductivity was only measured at the last round (+540 days), further interpretation of the resistivity changes regarding porosity, cementation, or clay content would be uncertain.Overall, this qualitative analysis of the monitoring results indicates that (a) the integral chargeability strongly increases at +40 days in both wells and tends to decrease at the following monitoring rounds, going back or close to the baseline, (b) the integral chargeability increase at +40 days is rather uniform over the whole depth interval in NN4 (200-400 m), while the increase in NN3 is localized between 450 and 550 m and (c) the resistivity tends to increase monotonously over time at the bottom (below 450 m) of NN3, which may correspond to a decrease of porosity related to intense pyrite and secondary mineral precipitation, although impossible to quantify further.
In the following section, spectral information is extracted from TDIP monitoring results to allow a first-order quantification of pyrite precipitation in areas with an increase in M int .

Extraction of Spectral Information Through Debye Decomposition
Since Debye decomposition tends to amplify the noise in the polarizability curves, decays that could not be fitted by a stretched exponential function (Equation 15) with R 2 ≥ 0.998 were discarded.This also smoothed out the decays that were kept and avoided fitting noise with the Debye decomposition.The theoretical justification and illustration of using a stretched exponential model for TDIP data is further described in Appendix B (Figures B1 and B2). Figure 4 also illustrates that the stretched exponential fit has minimal difference with the polarizability data in most cases.
Debye decompositions, illustrated in Figure 8 for two sets of decays in NN3 and NN4 (64″ spacing), show a similar trend as observed in Figures 6 and 7.The polarizability decays at +40 days are above the baseline as well as the later monitoring rounds, both in NN3 and NN4.This trend is also reflected in the coefficients m k , especially at long τ k values, and the resulting m tot , larger at +40 days.All Debye decompositions are available as Figures S1 to S4 in Supporting Information S1.
As can be seen on Figure B1, only 70% of the decays passed all the processing steps.Most affected data sets were at +380 days in NN4 (above 330 m) and at +40 days in NN3 (above 450 m).The noisy data above 450 m at +40 days is attributed to two factors that can be observed in Figure 3: (a) significantly stronger 50 Hz noise above 450 m than below at all dates (reason unknown) and (b) stronger 50 Hz noise at +40 days compared to baseline, at all depths.As can be observed in Figure 3, the signal level is always better with the 16″ (blue) than with the 64″ (red) electrode spacing, which is expected considering that the measured voltage decreases for increasing spacing with this configuration.This explains that fewer 16″ data were discarded during the Debye decomposition filtering procedure (B1).However, the 64″ data involve a larger depth of penetration around the borehole and are thus more relevant for conversion to pyrite volume change.
On the one hand, the integral chargeability presents the advantage of smoothing out noise and providing a qualitative overview of the trends for the whole data set.On the other hand, Debye decomposition allows extracting spectral information and further quantifying pyrite volume fraction change.Since the integral chargeability is the most easily measured parameter in the field, it would be practical to have a relation between m tot values obtained by Debye decomposition and the integral chargeability M int , as it would allow converting M int into pyrite volume fraction.However, we find that different relations exist between m tot and M int , depending on the mean relaxation time.Investigating further these relations is beyond the scope of this study.Thus, pyrite estimations in the next section are calculated on a reduced portion of the data set, where Debye decomposition was carried out.

Quantification of Pyrite Volume Fraction Change
The absolute difference in total chargeability, Δm tot , further converted into Δp v (Equation 14) are presented in Figure 9 for the three monitoring rounds.While the same trends are observed as in Figure 6, this conversion allows estimating the extent of pyrite precipitation or dissolution.Error bars also differentiate between significant differences in pyrite volume fraction versus differences too uncertain to be interpreted further.Overall, precipitation of up to 2% and 1% are observed with the 64″ and 16″ electrode spacing, respectively.Figure 9 illustrates that more volume fraction changes are observed when taking into consideration a larger volume around the boreholes, which emphasizes once more that changes are not limited to the very near vicinity of the injection boreholes.
In NN4, an increase in pyrite volume fraction at +40 days (+0.5%-1% compared to baseline) is observed with both the 16″ and 64″ electrode spacings at all depths, followed by a decrease at +380 and + 540 days.In addition, the 64" data indicate that the decrease at +540 days extends below baseline values ( 0.2%-0.5% compared to baseline).In NN3, changes are more depth-dependent than in NN4.Between 200 and 450 m, no significant changes are observed at any rounds.However, it is important to note that, at +40 days, most of the 64″ data in this depth-range were discarded during the Debye decomposition filtering procedure.Below 450 m, all the way down to 550 m, +1%-2% in the pyrite volume fraction is suggested by the 64″ data and +0.5%-1% by the 16″ data.At +380 and +540 days, the 64″ and 16″ data, where available, suggest a decrease toward baseline values, yet never extending below baseline values, as opposed to NN4.

Analysis of the Mean Relaxation Time
Overall, relaxation times are shorter in NN4 than in NN3, with a maximum of 300 msec in NN4 and above 1,000 msec in NN3 (Figure 10).It can also be noted that the variations of the average relaxation times across the different monitoring rounds follow the variations of m tot .In NN3 below 450 m, τ mean increases at +40 days and decreases back at +270 days.In NN4 across the whole well, τ mean increases at +40 days and decreases back at +270 days.In NN3, significantly longer relaxation times are observed below 450 m, which reflects the different  12); right: coefficients m k for all the τ k at the four measuring rounds.
shape of the polarizability decay curves, with most of the decay happening at late times, as opposed to depths above 450 m.Similarly, in NN4, decays at 346-348 m and at 380-382 m have a different shape and correspond to minimum and maximum values of τ mean over the whole well, 40 and 400 msec, respectively.
As was observed in Figure 1, longer relaxation times, corresponding to lower frequencies, can be associated with either larger metallic grain sizes or connected metallic grains, for example, forming veins.In this case, relatively small pyrite grain sizes are expected due to the short time available, especially after +40 days (around 10 μm according to laboratory experiments presented by (Prikryl et al., 2018)).Therefore, long relaxation times in NN3, below 450 m, may indicate connected pyrite particles in fractures.On the other hand, shorter relaxation times in NN4 may be associated with small, disseminated pyrite particles.While larger relaxation times may also be caused by larger fluid resistivity (Feng et al., 2020;Martin et al., 2022), very little resistivity changes are observed in NN3 and NN4 after +40 days so that the relaxation times should not be affected by fluid resistivity.Moreover, the larger relaxation times at the bottom of NN3 are associated to lower resistivity and therefore can only be attributed to the larger chargeability and grain sizes.
It is important to note here that TDIP data are only capturing polarization processes associated with relaxation times larger than 3 msec, which is above the typical range of relaxation times investigated by FDIP laboratory studies, such as Martin and Weller (2023) and (Revil, Abdel Aal, et al., 2015).Therefore, converting the relaxation time distributions into grain size distribution would be uncertain.

Interpreting Subsurface Processes Behind the Monitored Resistivity and TDIP Responses
In the field acquisition deployed here, a bulk TDIP response is measured, in which the various contributions to the chargeability (e.g., electrode vs. membrane polarization) are indistinguishable.As opposed to laboratory investigations in well-controlled experiments, conditions such as smectite content, porosity and fluid flow patterns are unconstrained and the interpretations thus rely on a series of assumptions.In particular, we assume that pyrite precipitation is the process governing the increase in chargeability.
The justification for converting total chargeability into pyrite volume fraction changes is discussed here.Based on reaction path modeling carried out in similar conditions, it is likely that some amount of smectite precipitates (Galeczka et al., 2022).Still, smectite precipitation is expected to have little influence on the IP signal in the presence of pyrite and could not alone explain the observed chargeability increases (Lévy, Gibert, et al., 2019;Lévy, Weller, & Gibert, 2019;Revil et al., 2017).Precipitation of iron-(oxyhydro)oxides following basalt dissolution is also possible, but only magnetic iron-oxides, such as magnetite or ilmenite, yield a strong polarization response (Peshtani et al., 2022).In these conditions, only non-magnetic Fe(III)-bearing phases are expected, such as goethite or hematite, based on laboratory experiments (Andreani et al., 2009;Menefee et al., 2018), field-scale modeling (Galeczka et al., 2022), and the general low-temperature alteration sequence of Icelandic basalt (Arnórsson et al., 1983;Crovisier et al., 1992;Gunnlaugsson & Arnòrsson, 1982;Stefánsson & Gíslason, 2001).In addition, amorphous silica and other silica-based minerals, such as quartz and zeolites, are  expected (Daval et al., 2011;Galeczka et al., 2022).Amorphous silica, zeolites, and Fe(III)-bearing phases further have a strong tendency to clog the pore space (Andreani et al., 2009;Menefee et al., 2018), without creating a distinctive polarization signal (Revil et al., 2002).Finally, the influence of microbial activity on the IP response is still poorly understood but seems to yield polarization of much smaller amplitude than the observed changes (Mellage et al., 2018) and will therefore not be considered here.
With these considerations in mind, TDIP monitoring brings three types of insights: spatial variability, long-term evolution, and quantity of pyrite volume changes.
The spatial variability of pyrite precipitation is assessed as a function of depth, variability between NN3 and NN4, and distance from injection well.First, it is clear that more pyrite precipitation occurs in the depth range 450-550 m in NN3 than in the rest of the well, while in NN4, rather uniform pyrite precipitation is inferred (Figures 6  and 9).The spatial variation in NN3 can be understood in light of the baseline resistivity (Figure 5), showing resistivity lower than 30 Ωm at levels deeper than 450 m in NN3.In this 450-550 m interval, more clay content was found at the time of drilling (Helgadóttir, 2021), indicating preferential hydrothermal flow.While the injection rate is ten times larger in NN4 than in NN3 (Table 1), more precipitation seems to occur at the bottom of NN3 compared to NN4.This may be due to fracture flow transporting the H 2 S-rich fluid rapidly away in NN4, not allowing time for H 2 S mineralization.There may exist an optimal permeability that maximizes the volumetric change in pyrite detected near the borehole, ensuring both ample fluid supply (H 2 S supply) and sufficient fluid residence time (basalt dissolution).
Second, it can be inferred that more precipitation occurs away from the well than on the borehole walls, based on the comparison between the 64″ and 16″ data sets in both wells (Figure 6).
The time-evolution of resistivity, integral chargeability, total chargeability and mean relaxation time show that subsurface processes may be more complex than the expected continuous pyrite precipitation suggested by previous studies (Clark et al., 2020;Gunnarsson et al., 2018).Integral chargeability time-series presented in Figure 7 show that in both wells, at most depths, the initial increase at +40 days is followed by a decrease, extending below baseline values in NN4 at depths below 330 m.This trend, which is common to the 6″" and 16″ spacings (Figure 6), is also observed with the total chargeability and relaxation time (Figures 10 and 9 and B1).
Here, three different interpretations are explored: (a) pyrite is still precipitating after +40 days, but the IP signal decreases anyway, (b) previously formed pyrite grains are "passivated" by a coating layer of secondary minerals preventing polarization while no more pyrite is forming and (c) pyrite is being re-dissolved.
Interpretation (i) would be consistent with the observation by Martin and Weller (2023) that the same amount of pyrite, with increasing grain size, yields a lower IP signal.However, the relaxation time, which is the main indicator for particle size (Figure 1) also starts decreasing after +40 days (Figure 10), which indicates that pyrite particles, if they are still there, are not growing in size.This interpretation is, therefore, ruled out.
Interpretation (ii) is based on the fact that strong polarization effects are caused by metallic particles only when they have a clean surface.If secondary minerals isolate the pyrite surface from the fluid, the particles are passivated, and the polarization decreases (Gurin et al., 2019).Amorphous silica and Fe(III)-bearing minerals are common coating layers, which may coat pyrite grain surfaces, as well as reduce basalt dissolution rates by coating fresh basalt and clogging the pore space (Andreani et al., 2009;Daval et al., 2011;Menefee et al., 2018).Corrosion of pyrite surface in the presence of little amounts of oxygen may also lead to a similar passivation effect, without significantly dissolving the pyrite particle (Placencia-Gómez et al., 2013).The monotonous resistivity increase at the bottom of NN3 tends to indicate progressive pore space clogging, which a series of nonpolarizable minerals may be responsible for.
Interpretation (iii) contradicts the common idea that pyrite mineralization is stable in this context once formed.However, practical field conditions differ from typical hypotheses in reactive transport models, which are constrained by the availability of suitable frameworks and databases.First, a high disequilibrium near the borehole, where H 2 S-rich fluids are constantly injected, makes conventional "steady-state" assumptions irrelevant.Second, previous experiments on CO 2 injection in basalt at Hellisheiði (similar set-up as the Nesjavellir field case presented here) triggered a bloom of iron-oxidizing bacteria, such as Gallionellaceae, and more generally sulfur-and iron-oxidation markers (Daval, 2018;Trias et al., 2017).At Nesjavellir, provided that oxygen was present in the system at some point (which could be the case due to e.g., leakage in the pipe system happening at the beginning of the H 2 S injection), specific microbial communities may have developed and triggered pyrite oxidation (Mielke et al., 2003;Percak-Dennett et al., 2017).The fracture networks that are ubiquitous in basaltic rocks could also be a pathway for intermittent oxygen delivery, which would extend the depth of the habitable zone for iron-oxidizing bacteria (Bochet et al., 2020).Overall, biologically-driven mechanisms, which could explain an unexpected pyrite dissolution over time, are not taken into account in previous studies claiming that pyrite mineralization is stable due to the difficulty of representing them correctly in reactive transport simulations.While high-temperature geothermal systems may prevent the development of most microbial communities, the shallow system studied here has a relatively low reservoir temperature, although the injected fluid temperature is higher (Table 1).
The pyrite volume fraction change estimated in Figure 9 assumes, to simplify, that all the decrease of chargeability is explained by pyrite dissolution.In practice, there may be a mix of passivation and re-dissolution and, thus, a more conservative interpretation of the vol.% presented in Figure 9 would be the change in pyrite volume fraction contributing to the IP signal.While it is impossible from geophysical data alone to determine whether passivation (interpretation ii) or re-dissolution (interpretation iii) are responsible for the observed signals, both could be further investigated by laboratory studies, as well as by installing a more continuous geophysical monitoring infrastructure.
In addition to the uncertainty on the source of polarization decrease, the quantification of pyrite volume presents additional uncertainty.It should first be noted that the influence of background polarization not related to sulfide minerals in the baseline (Figure 5) is managed by only converting the difference of total chargeability between different monitoring rounds into the difference in pyrite volume fraction.However, Equation 5proposed by Revil, Florsch, and Mao (2015) and used here presents other limits, as Martin and Weller (2023) found that the exact relation between m tot and p v depend on the grain size of pyrite particles.Furthermore, the use of an equation calibrated with FDIP measurements to convert TDIP measurements is made possible by "long" current injections (2 s), where the TDIP and FDIP responses are supposed to overlap (Martin et al., 2021).However, 2-s current injections remain short, compared to 64-s used by Martin et al. (2021) to show the equivalence with laboratory measurements.Therefore, the pyrite difference calculated here from 2-s TDIP data is certainly underestimated.It was also assumed here that, at a given depth, precipitation occurs uniformly over a certain (undetermined) lateral extent, considered representative of near-borehole conditions.Further modeling of the contribution from the different cylindrical layers around the borehole to the measured signal, and thus of the pyrite volume in each of these layers, could be achieved with inversion.
Finally, it is interesting to note that the mean relaxation time is not only time-dependent but also depth-dependent and that an overall consistency of the decay shapes at a given depth is seen over time (Figure 10).This suggests that the local lithology and fracture network influence the relaxation time, most likely by influencing the patterns for pyrite precipitation.Although this is purely speculative, it may have implications for a better understanding of fracture networks in subsurface storage systems.

Requirements and Added-Value of IP-Logging Compared to Current Practices
As mentioned in the Introduction, current monitoring of field-scale H 2 S reinjection, and similarly of CO 2 reinjection, have major pitfalls: (a) tracer tests capture less than 5% of the gas migration flow (Matter et al., 2016), (b) mass balance calculations assume that the amount of sulfur not captured by downstream measurements, compared to expectations based on tracer tests, has mineralized into sulfide minerals (Gunnarsson et al., 2018) and (c) longterm predictions by reactive transport models are not calibrated by field monitoring and thus present significant uncertainty (Aradóttir et al., 2012;Ratouis et al., 2022;White et al., 2020).Little has been researched on the longterm stability of the newly formed minerals in a context where the continuous injection of H 2 S (and similarly, CO 2 ) creates a constant disequilibrium in the system, in terms of acidity and redox conditions in particular.Overall, current monitoring approaches are still indirect and cannot physically confirm the presence, location, and quantity of expected minerals in the rock matrix, nor their long-term stability.The results presented here show the feasibility of monitoring pyrite precipitation with logging TDIP chargeability measurements within injection wells.However, several conditions/constraints/obstacles need to be overcome.
To start with, an appropriate tool is required.Compared to more conventional applications of TDIP, typically in mining exploration, the expected polarization signal is weaker (smaller volumes of sulfide minerals).At the same time, the presence of power lines due to the proximity to the power plant creates a significant 50 Hz background noise.Therefore, large current injection is essential to ensure a sufficient signal-to-noise ratio, and the acquisition of voltage with a high sampling rate is essential to separate the contribution from 50 Hz noise and from the rock to the voltage waveform.In addition, the need for quantification of pyrite precipitation, and thus comparison to frequency-domain IP calibrations, requires that current injection and voltage acquisition happen over a sufficiently long time to capture the polarization information needed for the interpretation.All these specifications required the development of new hardware in the existing QL40-IP instrument.In particular, the need for both longer injection/acquisition time and a higher sampling rate means that more data needs to be collected and recorded in the instrument's memory.The instrument went from recording 100 samples per cycle (depth) to 980.Nevertheless, a trade-off was adopted, and the current waveform, deemed less critical for the analysis, was not recorded; only the average current at each depth was recorded.In addition, the maximum injection/acquisition times were 2 s.Ideally, recording the current full waveform, injecting during longer times, and combining different injection times would allow a more accurate comparison to FDIP calibrations.
Field conditions are also important to consider to carry out successful TDIP logging.The instrument has a temperature limit of 70°C, and its exposure to corrosive fluids shall be limited.This can be challenging in the context of CO 2 and H 2 S injection wells at geothermal sites.On the other hand, injection wells are the location with strong changes in mineralogy.At the Nesjavellir site, the need to stop the injection and cool down the wells a few days before the measurements represented the main obstacle to repeating the measurements more regularly.In addition, the method requires an operator on-site to operate the instrument.In this context, where data redundancy is limited while TDIP measurements are easily subject to noise, data loss is a critical aspect to manage, especially for quantitative conversion.
Despite field limitations and remaining uncertainties, IP-logging monitoring brings novel insights into the processes at play, as emphasized in the previous subsection.A possible way forward would be to carry out the logging in dedicated monitoring wells relatively close to the injection wells to avoid the financial and environmental cost of stopping the injection during monitoring.Another possibility that would also fulfill the need for more continuous data acquisition to overcome episodic data losses would be the installation of permanent electrodes in dedicated wells.The electrodes would connect to an instrument remaining at the surface during the whole monitoring period, and only batteries would need to be changed periodically.Furthermore, this approach could allow 2D or 3D monitoring if several wells are used together in a cross-borehole electrical tomography approach.External data acquisition would also be helpful to further elucidate the physical processes behind the observed IP response.Monitoring the fluid conductivity over time appears essential to assess changes in fluid composition and further interpret resistivity changes.Complementary logging methods may include sonic and density logging for assessment of porosity changes, acoustic televiewer, and spinner logs for mapping fractures and flow paths, and color optical camera for direct observations of the precipitation of pyrite and biofilms at the borehole wall.

Can the IP-Logging Method Be Applied for Monitoring Other Reactive Processes?
Due to the strong signal caused by sulfides and other metallic minerals (Pelton et al., 1978), mapping ore deposits is still the main industrial application of IP methods today.Here, we present the first field example where reactive processes involving the precipitation and/or dissolution of metallic minerals are monitored with IP-logging.Other field examples for monitoring, with IP methods, reactive processes involving metallic minerals include the immobilization of uranium contamination by stimulating iron and sulfate-reducing microorganisms (Flores Orozco et al., 2011) as well as the injection of zero-valent-iron amendment to aid bioremediation of chlorinated solvents (Flores Orozco et al., 2015).In these examples, surface TDIP is used, and the limited spatial and temporal resolution complicates the interpretation.There, IP-logging could help improve the understanding of physical processes at play.
Studies on the link between IP and bacterial activity often find that microbially-mediated metallic minerals cause an IP response (Flores Orozco et al., 2011;Slater et al., 2007;Williams et al., 2005).However, recent studies suggest that polarization of bacterial cells themselves may control the IP response in certain cases (Mellage et al., 2018;Strobel et al., 2023).Therefore, in addition to monitoring the bioremediation of contaminants, IPlogging could also be applied to the study of microbial life in a dynamic natural environment, such as intermittent oxic-anoxic fluid mixing in fractured rocks (Bochet et al., 2020), provided that upscaling issues can be overcome (Mellage et al., 2018).
Finally, we evaluate to what extent the method developed here could be applied to monitoring CO 2 injection in basalts.Laboratory studies suggest that calcite precipitation causes significant polarization increase (Wu et al., 2010).However, others observe limited polarization increase, or even a decrease, in response to calcite precipitation (Wu et al., 2011;C. Zhang et al., 2012;Saneiyan et al., 2018) as well as no decrease associated with calcite dissolution (Halisch et al., 2018).Recent micro-and milli-fluidic investigations shed light on possible underlying IP processes at play.During calcite dissolution, the growth and transport of CO 2 gas bubbles create two-phase flow conditions, which may be responsible for a decrease of the quadrature conductivity, and a sequence of increase-decrease-increase of the in-phase conductivity (Rembert et al., 2023).During calcite precipitation, a strong increase in quadrature conductivity, only observed at the late stages of calcite precipitation, is attributed to pore space closing rather than mineral growth itself (Izumoto et al., 2022).The current petrophysical knowledge on the complex electrical response of calcite is insufficient to allow direct mapping of calcite in the field with IP methods.Nevertheless, H 2 S is often associated with CO 2 in emissions from geothermal or coal power plants (Fridriksson et al., 2016;Wang et al., 2011).H 2 S is often considered an impurity and is only injected with CO 2 to minimize the cost of separation.Still, the formation of pyrite could be used as a marker of mineralization processes, given that pyrite and calcite precipitations are primarily controlled by basalt dissolution rates (Prikryl et al., 2018;Stefánsson et al., 2011) and also share similar kinetics (Plummer et al., 1979;Williamson & Rimstidt, 1994).

Conclusions
To date, monitoring efforts of an implemented field-scale H 2 S reinjection system have only occurred via concentration measurements, mass balance calculations, and analysis of precipitates on submersible pumps in injection and downgradient boreholes.These studies, along with geochemical numerical models and laboratory simulations, concluded effective H 2 S sequestration through pyrite formation.Borehole sampling of fluids, although effective at identifying the occurrence of sulfide mineralization, provides limited quantitative, spatial, and temporal information on in-situ sulfide mineralization.
In our study, a "shallow" (200-550 m) injection of H 2 S at Nesjavellir geothermal field was monitored by means of resistivity and time-domain induced polarization wireline logging.Pyrite precipitation was expected to affect primarily the induced polarization signal, while other changes (e.g., fluid conductivity, temperature, clay minerals precipitation) were expected to affect the resistivity.Our results represent the first field demonstration of an IP monitoring in boreholes, designed to track mineralogical changes accompanying H 2 S injection into basaltic rocks, down to 550 m depth.
The QL40-IP logging tool, developed by the company Advanced Logic Technology, was modified to allow longer current injection and a larger sampling rate than previously.This was needed to filter a strong 50 Hz background signal, due to the proximity of a buried power line.The conversion of TDIP data into sulfide volumes was carried out by (a) developing a new strategy for outlier removal and smoothening IP decays, (b) applying an existing Debye decomposition approach to calculate the total chargeability and (c) convert the chargeability difference to the baseline into a volumetric change of pyrite.
The first logging measurements were carried out in the summer of 2019, and the baseline with the modified tool was measured in September 2020.H 2 S injection began on 29 January 2021.After the start of injection, five logging monitoring rounds (+40, +270, +380, +540, +630 days) were carried out in wells NN3 and NN4.Out of these five rounds, only three remained for each well due to too low current injected in the other rounds for unknown reasons.
Monitoring results show that polarization strongly increases at +40 days in both wells and tends to decrease at the following monitoring rounds, going back close to the baseline or below.The polarization increase is rather uniform over the whole depth interval in NN4 (200-400 m), while the increase is clearly localized between 450 and 550 m in NN3, where more hydrothermal alteration was also found prior to the injection of any wastewater.Conversion of chargeability absolute difference into pyrite volume fraction change indicates precipitation of up to 1% in NN4 and 2% in NN3 at +40 days.In both wells, changes are more pronounced with the larger electrode spacing (64"), indicating that pyrite precipitation takes place away from the well.
The decrease of chargeability observed in both wells at +270 and +380 days suggests that pyrite is either passivated or re-dissolved after precipitating.Overall, TDIP monitoring is sensitive to a sequence of overlapping reactive processes resulting from H 2 S injection, showing its complementarity to traditional geochemical montioring.Recommendations for future monitoring include: (a) installation of permanent monitoring systems to better understand the trends in an in-situ controlled experiment, (b) considering the possibility of pyrite redissolution in reactive transport modeling of H 2 S sequestration, and (c) validation of geophysical results by additional in-situ measurements and observation, such as color camera and multi-level water sampling in wells, as well as by post-injection core retrieval for physical evidence.The presence of specific microbial communities that can produce oxygen at these depths should also be further investigated.
Overall, our results provide pieces of answers and insights into the subsurface physical processes controlling mineralization and its sustainability over time, as well as quantitative information on the sulfide mineralization with high vertical resolution.While we were not able to bring definitive answers on the stability of pyrite over time, we identified a previously unknown potential problem to H 2 S storage approaches.Future research exploring to what extent the IP response of re-dissolution and passivation may be discriminated would be valuable.not analytical since the HN and the KWW relaxation functions are not exactly Fourier transforms of each other.In particular, the relation between the relaxation times in both functions depends on the exponent beta in the KWW function, β.For example, for β close to 1, τ HN ≃ τ KWW , and for β = 0.1, τ HN = 300τ KWW .The KWW has one less parameter, so it can always be transformed into an HN model, but the contrary is not always true.The product α HN γ HN is close to β.In addition, even though the HN model is a generalized version of the Pelton model, the fitted chargeability parameter cannot be used to quantify pyrite as it is not strictly comparable to the chargeability in the Pelton model.
That's why the Debye decomposition is used here to relate the measured TDIP data with petrophysical relationships calibrated on FDIP data.
Nevertheless, based on the theory described by Alvarez et al. (1991), the stretched exponential function is used to automatically discriminate noisy discharge curves that do not follow the stretched exponential shape.Debye decomposition is only applied to discharge curves that can be reasonably fitted by a stretched exponential.Moreover, the stretched exponential fit allows smoothing out the "accepted" discharge curves that may still contain some noise, which in turn results in more meaningful outcomes from the Debye decomposition.This is illustrated for the total chargeability, M tot , and the mean relaxation time, τ mean in Figures B1 and B2.In these figures, Debye decomposition carried out on the stretched exponential fitted function is compared to Debye decomposition carried out on the polarizability data themselves (before the fitting step).The comparison clearly indicates less noise in the curves of SE fitted spectra.

Acknowledgments
The logging crew at ÍSOR is thanked for the numerous logging campaigns: Halldór Ö.

Figure 1 .
Figure 1.Laboratory results for naturally altered volcanic samples from the Krafla geothermal field (north-east Iceland) containing varying amounts of pyrite with different grain sizes.FDIP data for four different samples: (a) impedance modulus and (b) phase angle, as a function of frequency.(c) Scanning electron microscope (SEM) images for the three samples containing pyrite.Pyrite grains are highlighted in red.(d) Distribution of pyrite grain size, based on the SEM images of the three samples.
Two injection wells are used for monitoring: NN3 and NN4, drilled into fresh lavas to 563 m (NN3) and 422 m (NN4) depth in 2001, northeast of the power station.Their specific lithology and alteration stage at the time of drilling (i.e., prior to large-scale warm wastewater injection) are presented inHelgadóttir (2021), emphasizing an increase in the amount of clay minerals below 400 m in NN3.Between 2004 and 2021, injection of the initial wastewater into the cold groundwater system took place through these wells and has probably been enhancing basaltic rock alteration compared to natural conditions.Starting from January 2021, H 2 S has been continuously injected through the final wastewater.

Figure 2 .
Figure 2. Aerial view of Nesjavellir valley (south-west Iceland).The white building shows the geothermal power plant.The two wells, NN3 and NN4, investigated in this study are indicated by white circles, while other wells in the area are shown with small black rectangle signs.A buried power line (132 kV), most likely responsible for the 50 Hz background noise, is shown in blue.Pipelines are shown in pink and roads in black.

Figure 3 .
Figure 3. (a) Full waveform (FW) TDIP data measured in NN3 with 16″ and 64″ electrode spacings.The 50 Hz noise is further illustrated in (b) by zooming on the positive decay for the 64″ data (signal in black), where the FW signal after removing the 50 Hz noise is also shown (in red).Further processing of the FW data into polarizability curves (after dividing by DC voltage) is shown in (c) with linear scale and (d) with log-log scale.Results are shown at two different dates (baseline and +40 days) and two different depths (263 and 470 m), showing varying levels of background noise.

Figure 4 .
Figure 4. Resulting polarizability decays after the five processing steps in (a) NN3 and (b) NN4.An example of noisy decay in NN3 is shown at 382-384 m depth, with larger error bars and a non-exponentially decaying shape.

Figure 5 .
Figure 5. Resistivity (blue/red) and integral chargeability (green/black) logs in (a) NN3 and (b) NN4 using the 16-and 64" spacing (colors).All plots are based on an averaged signal every 2 m, and the error bars correspond to the standard deviation within these 2-m intervals.The dashed line shows the resistivity threshold of 30 Ωm.The gray rectangles indicate intervals with locally high resistivity and low integral chargeability.

Figure 6 .
Figure 6.Relative difference of the (a) and (c) temperature-corrected resistivity ρ and (b) and (d) integral chargeability M int in NN3 (top) and NN4 (bottom), for both 64" and 16" electrode spacing.Colors indicate the monitoring rounds.

Figure 7 .
Figure 7. Evolution of the integral chargeability (absolute difference in mV/V) over time in (a) NN3 and (b) NN4 in different depth intervals, represented by the color scale.The dashed black line indicates the baseline level.Note the different scales for NN3 and NN4.Only data measured with 64″ spacing are shown.

Figure 8 .
Figure 8. Debye decomposition for two sets of decays measured with the 64″ spacing in (a) NN3 at 492-494 m and (b) NN4 at 362-364 m.Left: polarizability data and fit by the RTD function (Equation12); right: coefficients m k for all the τ k at the four measuring rounds.

Figure 9 .
Figure 9. Estimations of pyrite fraction change (absolute difference of pyrite vol-%) in (a) NN3 and (b) NN4, based on TDIP logging with 64" and 16" electrode spacing, using Debye decomposition of the decays.

Figure B1 .
Figure B1.Total chargeability from Debye decomposition in NN3 (bottom) and NN4 (top), with 64″ (left) and 16″ (right) spacing.The results are shown for Debye decompositions carried out on stretched exponential fit ("SE fit") and on original data.Note that the Debye decomposition on original data includes more discharge curves since it happens before the filtering procedure related to the SE fit.

Figure C1 .
Figure C1.Temperature measured in (a) NN3 and (b) NN4 at the four monitoring rounds.
Starting in 2004 and prior to 2021, wastewater composed of condensate and separation water has been disposed of in drainage wells.This initial wastewater had negligible concentrations of dissolved CO 2 and H 2 S. Starting on 29 January 2021, a third type of fluid, called seal water, is now mixed into the condensate and separation water before injection.Seal water contains the non-condensable fraction of gases, primarily H 2 S, that pose corrosion problems to turbines over time.It is formed by dissolving the CO 2 and H 2 S gases from liquid ring vacuum pumps at the condensers in cold groundwater or condensate water.The currently injected wastewater contains, on average, 75 ppm of H 2 S (Table1), with the intention to sequester H 2 S through mineral storage.

Table 1
Details on the Nesjavellir NN Injection Wells in This Study Note.The borehole fluid conductivity was measured with a QL40-FTC logging probe (Mount Sopris QL40-FTC, 2014) in July 2022 (+540 days).Temperature and caliper logs are presented in Figures C1 and C2 in Appendix C.

Table 2
TDIP Monitoring Timeline and Evaluation of Data Quality LÉVY ET AL.
Stefánsson, Halldór Ingólfsson, FriðgeirPétursson.Jakob J. Larsen is thanked for their very appreciated help with signal processing.Thus S. Bording, Pradip K. Maurya, and Gianluca Fiandaca are thanked for all the valuable discussions around the processing and interpretation of TDIP data.Ásdís Benediktsdóttir, Vala Hjörleifsdóttir, Íris E. Einarsdóttir and Bjarni S. Gunnarsson are thanked for their help in organizing and coordinating logging campaigns with the power plant.Iwona Galeckza, Helga M. Helgadóttir, and Sveinborg H. Gunnarsdóttir are thanked for numerous discussions around the geology and geochemistry at Nesjavellir.LL received funding for the GEMGAS project from the Technological Innovation Fund, Taekniþróunarsjóður, at the Icelandic Center for Research, Rannsóknamiðstöð Íslands(Grant 198637- 0611).DC received a PhD grant from the Nordic Center for Volcanological research, Nordvulk.AW acknowledges the support of Deutsche Forschungsgemeinschaft through Grant 425975038.