Water–rock interactions in the Bruchsal geothermal system by U–Th series radionuclides

Uranium and thorium decay series disequilibria in deep geothermal brines are a result of water–rock interaction processes. The migratory behavior of radionuclides provides valuable site-specific information and can therefore be an important tool for reservoir characterization and sustainable management of geothermal sites. In this study, we present data from long-term monitoring of naturally occurring 238U, 232Th and 235U series radionuclides analyzed in brine samples collected from the Permo-Triassic sedimentary reservoir rock at the Bruchsal geothermal site (SW Germany). The results show that radionuclides of the elements radium (226Ra, 228Ra, 224Ra, 223Ra), radon (222Rn), and lead (210Pb, 212Pb) are rather soluble in brine, while isotopes of uranium (238U, 234U, 235U), thorium (232Th, 228Th, 230Th), polonium (210Po), and actinium (227Ac, 228Ac) have low solubilities and are mostly immobile. Activities of radium isotopes in the geothermal brine exceed those of their thorium progenitors (average 226Ra = 29.9 Bq kg−1, about 103 times that of its 230Th parent). Modelling the observed disequilibria allows the following conclusion on water–rock interaction processes: (1) supply from alpha-recoil depends on isotope half-life because it is limited by the rate of diffusion through microfractures causing isotopic fractionation. (2) Radium retardation due to adsorption is low (226Ra/222Rn = 1.3) resulting in adsorption–desorption rate constants in the order of 10−10 s−1 for k1 and 10−9 for k2. (3) Scavenging of 226Ra from brine can best be explained by co-precipitation with barite resulting in an observed 226Ra anomaly in the solids of the reservoir section. The precipitation rate constant amounts to ca. 3.4 × 10−8 s−1 corresponding to a mean removal time of radium from brine by mineral precipitation to approximately 1 year.

brines interact with solid phases with whom they come into contact. In consequence of these water-rock interactions, an elemental fractionation occurs resulting in a state of disequilibria (Osmond and Cowart 1992). Such radioactive disequilibria were found in deep geothermal brines in the Upper Rhine Graben (URG). Here, radium isotopes ( 228 Ra, 226 Ra, 224 Ra, 223 Ra) have activities far exceeding those of their thorium progenitors (Eggeling et al. 2013). Radioactive disequilibria caused by the preferred solution of radium were also documented for deep geothermal brines in the Salton Sea Geothermal field (Zukin et al. 1987). Previous studies have shown that radium concentrations are often high in saline waters (Kraemer and Reid 1984;Dickson 1985) and geothermal brines (Hammond et al. 1988;Rihs and Condomines 2002;Condomines et al. 2012), but rather low in low-temperature, low-salinity groundwaters (Krishnaswami et al. 1982;Luo et al. 2000;Porcelli 2008).
Modelling of these disequilibria provides information about the respective water-rock interaction processes controlling radionuclide supply into and scavenging from solution, respectively. This information is very useful to investigate the long-term migratory behavior of uranium and thorium series radionuclides which is not only an important issue for the integrity of nuclear waste disposals, but also for geothermal sites in terms of handling radionuclide-bearing precipitated minerals (scales) in surface installations: natural radionuclides which are once mobilized by water-rock interaction processes in the reservoir may be transported with the fluid through the geothermal power plant and trapped in solid solutions because of changing temperature and pressure conditions. Furthermore, understanding the migratory behavior of radionuclides in the reservoir may be useful for the characterization and sustainable management of the geothermal reservoir. Tricca et al. (2001) describe the water-rock interactions as physico-chemical reactions between three phases: the aqueous phase, the solid minerals and a reactive surface layer with specific properties, area and thickness. The transfer rate of a radionuclide from the rock material into solution depends on: (a) the in situ radioactive decay of its dissolved Fig. 1 Uranium-radium, thorium and actinium decay series. Vertical arrows indicate alpha decay, diagonal arrows indicate beta decay. For radionuclides with more than one decay mode, only the most frequently occurring is given parent; (b) the desorption from the surface coating; (c) the alpha-recoil across the solidliquid interface within a distance of several tens of nanometers, and (d) the dissolution of the aquifer solid. The removal of a radionuclide from the brine depends on: (a) its radioactive decay in solution; (b) the adsorption onto the surface layer, and (c) the incorporation into precipitates.
In the past, several mathematical solutions of simplified aquifer models dealing with naturally occurring radionuclides have been developed. Andrews et al. (1982Andrews et al. ( , 1989 considered physico-chemical mechanisms for radionuclide supply and removal and calculated the timescale of water-rock interactions. However, the authors did not consider transport by advection. Krishnaswami et al. (1982) computed rate constants of sorption processes. Furthermore, they determined residence times of daughter nuclides by means of alpha-recoil input from 222 Rn activities, although without considering the effects of advective transport and mineral dissolution/precipitation. Davidson and Dickson (1986) suggested a model of uranium and radium isotopes transport including dispersive flow, but without considering precipitation and dissolution processes. Ku et al. (1992) proposed a model that accounts for radionuclide transport by advection and first-order kinetics, sorption-desorption, dissolution-precipitation of U-Th series radionuclides processes as well as the supply from alpha-recoil. Tricca et al. (2000Tricca et al. ( , 2001 suggested a model for the combined groundwater transport of naturally occurring U, Th, Ra, and Rn isotopes with regard to advective transport as well as the physico-chemical processes of weathering, decay, alpha-recoil and sorption at the water-rock interface. The present study is aimed at better understanding of the behavior of uranium and thorium series radionuclides in the Bruchsal geothermal brine. A comprehensive dataset was generated by frequent, long-term fluid sampling. Results of the geochemical surveys of major/minor elements as well as isotopic measurements of U, Th, Ac, Ra, Rn, Po, Bi and Pb are presented below. From the observed isotope disequilibria, water-rock interactions were investigated and their effects on radionuclide transport in the geothermal reservoir assessed. Since the modelling of radioactive disequilibria requires both types of information, the composition of the fluid and the solid phase, the authors refer to their previous work (Kölbel et al. 2020) where the Bruchsal reservoir rock was intensively examined based on cuttings from the geothermal boreholes. The modelling part of this study focuses on radium isotopes since their range in half-lives and their interrelation in the respective decay chain (cf. Fig. 1) allows the determination of water-rock interaction rates across different timescales.

Upper Rhine Graben (URG)
The Upper Rhine Graben is part of the European Cenozoic Rift System that extends from the Mediterranean to the North Sea coast (Ziegler 1992). It is characterized by an NNE-SSW striking extension structure with a length of approximately 300 km and a width of up to ca. 40 km. The deep Hercynian basement consisting of Paleozoic granites is overlain by clastic sediments (sandstones) from Permian to Lower Triassic and by Middle Mesozoic to Cenozoic sediments. The base of the sediments in the center of the valley is ca. 3000 m deeper than at the valley shoulders (Ziegler 1992).
The Upper Rhine Graben offers favorable conditions for the exploitation of geothermal energy. This is supported by spatially varying local heat flow anomalies and temperature anomalies at large depths. Most of the thermal anomalies are related to large-scale fluid circulation (Pribnow and Schellschmidt 2000).

Bruchsal site (Germany)
The Bruchsal geothermal site is located at the eastern main boundary fault of the Upper Rhine Graben. The geothermal power plant consists of a borehole doublet: an injection well (GB1) and a production well (GB2) located at a distance of 1.5 km from each other. Because of the complex tectonic structure, the reservoir section at the injection well differs from that of the production well with respect to depth and thickness.
The geothermal reservoir is located in Permo-Triassic sedimentary rocks, characterized by large-scale normal faults of varying step heights, ranging between 20 and 350 m (Fig. 2). The hydrothermal reservoir is developed in a horizon comprising Middle Buntsandstein to Upper Rotliegend rocks (depth interval: 2220-2485 m). The main inflow zones are located in the fractured zones of the Upper Permian at the depth interval between 2440 and 2470 m (Joachim et al. 1987). Hydrothermal alteration processes resulted in an almost complete transformation of the original feldspars  Joachim et al. (1987) and Kölbel et al. (2020)

Sampling and analytical methods
Brine samples were collected at a sampling point close to the GB2 well head at the Bruchsal geothermal site. Hydrogeochemical investigations were conducted on thermal water abstracted from the production well GB2 during fluid circulation. The sampling period spanned between October 2016 and May 2017. In total 32 water samples were collected. Major ions and selected trace elements were analyzed in the samples by ICP-OES, ion chromatography and photometry. Temperature, pH and conductivity were measured on-site. 238 U, 232 Th and 235 U decay series radionuclides ( Fig. 1) were analyzed in 13 samples well GB2 to study their temporal variations in activities. The main research focus was on activity variations of radium isotopes ( 224 Ra, 223 Ra, 228 Ra, 226 Ra). Radium-sampling was conducted using gas-tight 1.2-l Marinelli beakers (type G-130 G) prepared with 2.5 ml 65% HNO 3 (suprapur) to inhibit precipitation of solids. Water samples were not filtered.
Gamma spectrometry surveys were carried out using a p-type HPGe coaxial detector of 30% efficiency (with respect to 3″ × 3″ NaI(Tl) detector). The germanium crystal had a diameter of 76 mm. The detector was embedded in a 10-cm lead shield to protect against background radiation. Specific activities of radium isotopes were determined by gamma spectrometry allowing simultaneous measurements without further sample preparation. The list of the gamma rays used for the determination of activities of radium isotopes is reported in Table 1. Measurements (M t1 ) were performed immediately after sample collection to determine activities of the shortlived radium daughters ( 228 Ac, 212 Pb, 214 Pb). The measurement duration ( t ) was ca. 80,000 s (= 22.2 h). A second measurement (M t2 ) was carried out after storing the sample for more than 42 h, but less than 350 h.
In addition, radionuclides of the 238 U, 232 Th and 235 U decay series were analyzed in an external certified laboratory. The activities of 238 U, 234 U, 230 Th, 226 Ra, 210 Po, 235 U, 227 Ac, 223 Ra, 232 Th, 228 Th and 224 Ra in the Bruchsal brine were measured by alpha spectrometry. 228 Ra was measured by beta counting. The accuracy of the gamma spectrometric method was checked by alpha spectrometry. Moreover, activities of uranium and thorium isotopes were determined which is not possible by gamma spectrometric measurements alone due to low U-Th activities in the brine samples.
Reflecting the sample point at well head, the radionuclide activities measured in the lab were corrected by considering the travel time of produced fluids from reservoir to surface (lag time correction) as well as the elapsed time since fluid sampling.
Brine density was calculated according to Mao and Duan (2008), considering reservoir temperature and pressure as well as the molality of the NaCl brine (T = 135 °C, p = 250 bar, M(NaCl) = 2.1 mol/kg). The resulting brine density is 1023.24 kg/m 3 . The activities of dissolved radionuclides are reported as disintegration rate per fluid-mass (atoms s −1 kg −1 ).

Ra
The 226 Ra activity was directly determined using its gamma ray energy at 186.2 keV. A possible interference with 235 U (185.7 keV) is negligible due to the low uranium activity shown by high-resolution ICP-MS. Alternatively, the gamma-peaks of 214 Pb and 214 Bi (daughters of 226 Ra and 222 Rn, respectively) can be used to calculate the 226 Ra activity. In this case, the measurement can be performed after 20 days at the earliest, because this time is required to reach secular radioactive equilibrium (assuming no radon loss).

Ra
154.3 keV gamma ray energy was used to determine the 223 Ra activity. Here the presence of an interfering peak from 228 Ac has to be considered (cf. Table 1). According to Condomines et al. (2010), 223 Ra activity based on the 154.3 keV peak can be corrected by: where 223 Ra and 228 Ac are average activities integrated over the counting time. C and C S are the counts for the sample and the standard, respectively. ε Ac and ε Ra are the (1) 223 Ra = C C S · 223 Ra S + 228 Ac S · ε Ac I γ , Ac ε Ra I γ , Ra − 228 Ac · ε Ac I γ , Ac ε Ra I γ , Ra ,  Fig. 1).
Measurements of both Ra isotopes include the count rates of their short-lived daughters 228 Ac (t 1/2 = 6.13 h) and 212 Pb (t 1/2 = 10.64 h) which evolve through time during counting. The time-dependent evolution of the thorium decay series radionuclides is illustrated in Fig. 3 and can be described by a system of coupled differential equations. Their general solution was first given by Bateman (1910). A radioactive decay chain (N 1 → N 2 → … → N i ) with the decay constant λ i can be described by the following differential equations: Assuming zero concentrations of all daughters at time zero, Bateman (1910) expressed the concentration of nth radionuclide after time t as: where the coefficients are calculated by Degering and Köhler (2011) adjusted Eq. (4) to the time-averaged activity expressed as: where The mathematical solution for the Bruchsal site was computed by employing Mathcad ® . For better understanding, index notation of 228 Ra = A 1 ; 228 Ac = A 2 , 228 Th = A 3 and 224 Ra = A 4 are used from here on.
The activity of 228 Ra at sampling time (t = 0) was determined from the 228 Ac activity after a waiting period t d of minimum 42 h after sampling: 224 Ra was determined by gamma-rays emitted by the daughter nuclides 212 Pb (238.6 keV) and 208 Tl (583.1 keV) with the 212 Pb/ 224 Ra ratio reaching a steady-state value of 1.14 after ca. 100 h:

Major and trace elements
Field parameter and major and minor element data are reported in Table 2. The Bruchsal brine is highly concentrated in chloride, sodium and other alkali metals and alkaline earth metals, containing up to 131 g/l of total dissolved solids (TDS). Furthermore, the brine is enriched in sulfate and hydrogen carbonate as well as heavy metals such as lead, arsenic, and cadmium. In contrast, the concentration of organic compounds is low. Eh-pH conditions are difficult to determine because of the change in pressure and (5) temperature between reservoir and the sampling location at ground level. At the sampling point, pH values range between 5.0 and 5.9 [], while Eh values are relatively constant at ca. 81 mV (on average).
Activities of Th ( 230 Th, 232 Th and 228 Th) and U ( 238 U, 234 U, 235 U) isotopes in brine were below the limit of analytical determination. This suggests that in the reservoir uranium preferentially exists in the tetravalent state, forming insoluble phases such as UO 2 and USiO 4 . Thorium is only stable in the tetravalent state; irrespective of the redox conditions, uranium requires a reducing environment to become tetravalent as in deep geothermal reservoirs (Attendorn and Bowen 1997).
In contrast, high levels of activity were observed for Ra isotopes ( 226 Ra, 228 Ra, 224 Ra, 223 Ra), 222 Rn and Pb isotopes ( 210 Pb, 212 Pb). Their activities are several orders of magnitude higher than those of their thorium progenitors. 222 Rn, is in slight excess relative to its parent 226 Ra.
Isotopes of Pb ( 210 Pb and 212 Pb) have high solubilities resulting in fluid activities > 15 Bq kg −1 . 210 Pb (A 210Pb = 25.8 Bq kg −1 ) was found to be deficient relative to its progenitors 222 Rn (A 222Rn = 37.8 Bq kg −1 ) and 226 Ra (A 226Ra = 29.0 Bq kg −1 ), but still in  the same order of magnitude indicating that Ra, Rn and Pb have a comparable mobility in the Bruchsal geothermal system. The activity of 210 Po is below the limit of determination. The low 210 Po activity in comparison to its parent 210 Pb suggests removal of 210 Po from the brine. The same applies to isotopes of actinium: the short-lived 228 Ac (A 228Ac = 6.1 Bq kg −1 ) of the thorium decay series is deficient relative to its parent 228 Ra (A 228Ra = 15.9 Bq kg −1 ). 227 Ac, an isotope of the 235 U decay series, was found to have a lower activity than its daughter 223 Ra indicating the rapid depletion of 227 Ac from brine.
The rock/brine activity ratio (R c ) is a measure of the relative mobility of the isotopes (Zukin et al. 1987). Values for R c were determined based on analyses of brine (Table 3) and rock samples from borehole cuttings (Kölbel et al. 2020). The results are summarized in Table 4. High R c values in the order of magnitude of 10 2 for isotopes of U ( 238 U, 234 U, 235 U) and Th ( 232 Th, 228 Th, 230 Th) reflect the affinity of the nuclides to the surface of the solids for both elements. In turn, isotopes of Ra ( 226 Ra, 228 Ra, 224 Ra, 223 Ra) and Pb ( 210 Pb and 212 Pb) have lower R c values and therefore, a relatively high mobility in the investigated geothermal system is implied.
228 Th (a decay product of 232 Th) has a rather low R c value compared to 232 Th. This reflects a better accessibility of 228 Th to the geothermal brine because of the good solubility of its 228 Ra progenitor. Following the decay of 228 Ra, 228 Th is adsorbed onto grain surfaces producing a comparatively high activity compared to the activity of the 232 Th isotope.
In summary, isotopes of thorium, polonium and actinium generally display low levels of activity in the brine as a result of their poor solubility. The processes of removal of the above isotopes from solution are most likely adsorption and/or co-precipitation.

Table 4 R c values of Th-U decay series radionuclides
R c is defined as the ratio of rock activity relative to average fluid activity (cf.  Ra,224 Ra and 223 Ra activities show lower values that vary between 15 and 17 Bq kg −1 , 9 and 12 Bq kg −1 , and 0.3 and 0.7 Bq kg −1 , respectively. Variation in activity levels of the radium isotopes are likely a consequence of analytical uncertainties and steady state in activities of the radium isotopes can be assumed.
This observation is also made with the activity ratios of radium that are principally constant within 2σ analytical uncertainties during the period of sampling: 228 Ra/ 226 Ra, 224 Ra/ 228 Ra and 223 Ra/ 226 Ra ratios display mean values of 0.55 ± 0.07, 0.65 ± 0.07 and 0.02 ± 0.01.
Gamma spectrometry results from earlier analysis (cf . Table 6), however, indicate that Ra activity in brine has increased over the past years. This particularly applies to 226 Ra and 228 Ra whose activities both increased by 25% resulting in a constant 228 Ra/ 226 Ra ratio (0.53 ± 0.09 in 1986 and 0.55 ± 0.07 in 2016/17). However, the short-lived 224 Ra activity has been more or less stable since 1986 resulting in a decrease in the 224 Ra/ 228 Ra activity ratio over the past 30 years from 0.75 in 1986 to 0.65 in 2016.

Modelling the disequilibria in water-rock systems
Water-rock interaction processes in the Bruchsal geothermal reservoir were mathematically modeled based on Ku et al. (1992). This model allows for physico-chemical reactions as well as advective transport. The model is robust in terms of input parameters and focuses on the simulation of the behavior of radium in the rock-brine environment. Employment of the various radium isotopes provides a way of quantifying relevant parameters of the water-rock system due to their wide range of half-lives. Theoretically, the model can be applied to all elements with numerous instable isotopes such as U ( 238 U, 235 U, 234 U) and Th ( 234 Th, 232 Th, 231 Th, 230 Th, 228 Th, 227 Th). However, Th-U isotopes are not used in this study due to their poor solubility in the Bruchsal brine and the resulting lack of data.

Model assumptions suggested by Ku et al. (1992)
The processes of sorption-desorption and dissolution-precipitation of radionuclides are determined by reaction kinetics. Ku et al. (1992) subdivided three "pools" in which radionuclides can reside: the dissolved, adsorbed and solid pool. Figure 5 depicts a schematic representation of the conceptual model including water-rock interaction  processes and the three different "pools". Model parameters are listed in Table 7. Ku et al. (1992) defined the following model assumptions: (1) In the dissolved "pool", radionuclides are exchangeable with those in the adsorbed pool, but not with those in the solid "pool".
(2) Transfer of radionuclides between the dissolved and solid "pools" is achieved in particular by dissolution, co-precipitation and alpha-recoil. (3) Dissolution and precipitation are considered irreversible because dissolved nuclides have limited and very slow communication with the solid "pool" which is located further inside the rock matrix. (4) Alpha-recoil input from the adsorbed and dissolved "pools" to the solid "pool" is negligible. (5) Distributions of radionuclides in solid, adsorbed, and dissolved "pools" remain stationary.

Governing equations
Based on mass balance, the activity of a given radionuclide dissolved in a volume of brine with a constant density can be expressed as (Luo et al. 2000): where Q is the supply rate by water flow, atoms kg −1 s −1 ; P w is the supply rate of radionuclide to fluid by dissolution, atoms kg −1 s −1 ; P r is the supply rate of radionuclide to fluid by alpha-recoil, atoms kg −1 s −1 ; R f is the retardation factor due to adsorption and desorption, dimensionless; A is the radionuclide activity in brine, atoms kg −1 s −1 ( = C ); k p is the first-order precipitation rate constant, s −1 ; C is the radionuclide concentration in brine, atoms kg −1 ; and ′ is the superscript referring to the radioactive parent.
The retardation factor R f is formulated as follows (Krishnaswami et al. 1982): where is the radioactive decay constant of radionuclide, s −1 ; K is the dimensionless distribution coefficient; k 1 is the first-order adsorption rate constant, s −1 ; and k 2 is the first-order desorption rate constant, s −1 . For all radium isotopes whose thorium parents are quite insoluble in the geothermal fluid, R ′ f A ′ is negligible: Furthermore, processes of dissolution and precipitation will not influence the activity of short-lived radionuclides and thus, by setting P W = 0 and k p C = 0, Eq. (16) may be simplified to (17) Q + P r = R f A, C i and C are initial and measured concentrations, respectively, and τ b is the transit time of brine in the aquifer. Positive or negative values of Q denote net gain or loss due to fluid flow, i.e., advective transport (Luo et al. 2000).

Alpha-recoil
Alpha-recoil describes a process in which a radioactive daughter is mobilized from its initial position by the energy of an alpha decay (Sun and Semkow 1998). During the decay, an atomic nucleus emits an alpha particle. The released ionizing radiation has an energy content of 4-6 MeV. Because of the law of conservation of momentum, the emitted alpha particle and recoiling nucleus will each have a well-defined energy after the decay. Because of its smaller mass, most of the kinetic energy is transferred to the alpha particle. The recoiling nucleus will have a kinetic energy in the order of 100 keV (Sun and Semkow 1998). Nevertheless, the energy transfer to the decay product is high enough to shift atoms that are close to the mineral surface out of the mineral grain into the pore space (Fig. 6). The probability of recoiling out of mineral grains depends on the isotope recoil distance and the number of previous alpha decays. While 228 Ra is directly formed by the decay of 228 Th, 226 Ra is formed by three alpha decays of 238 U. Thus, the probability of 226 Ra to end up in the fluid by alpha-recoil is significantly greater than the probability for 228 Ra due to the greater amount of the three alpha energies.
Alpha-recoil supply rates ( P r ) for radium isotopes can be estimated from the activities of the decay series progenitors ( 238 U, 235 U and 232 Th) in the adjoining rock expressed as (e.g., Kigoshi 1971): where A ′ is the parent activity in solids (atoms kg −1 s −1 ), r is the recoil distance (Å), ε i is the recoil efficiency for nuclide i (), S is the surface area of solids, expressed as area of solid per volume of fluid contacting the solids (m 2 m −3 ) also called flow-wetted surface and ρ s is the density of solid (kg m −3 ).  (Kölbel et al. 2020). Surface area per mass is about 2000 m 2 kg −1 for the Permo-Triassic sandstones (Heap et al. 2019). This corresponds to a flow-wetted surface of 1.0 × 10 8 m 2 m −3 based on a porosity of 0.05. Sun and Semkow (1998) published data of Ra recoil distances Retardation factor due to precipitation as well as adsorption and desorption -  (Table 8).
Since recoil efficiencies for radium isotopes are not that easy to determine, ε i is part of the discussion.

Retardation factor and distribution coefficient
The retardation factor R f describes the flow rate of the fluid relative to the rate of migration of a radionuclide in the flow (Ku et al. 1992). The separation between the adsorbed and dissolved nuclides through chemical exchanges might be stated by the dimensionless distribution coefficient K (Krishnaswami et al. 1982): where A a is the activity of an adsorbed radionuclide (atoms per equivalent fluid-volume) and A is the radionuclide activity dissolved in solution (atoms per fluid-volume).
Since distribution coefficients are usually determined by adsorption-desorption experiments in the laboratory and therefore, expressed in units of volume per mass, K can be derived from where K d is the distribution coefficient, volume mass −1 ; ρ s is the density of aquifer solids, mass solid volume −1 ; and φ is the porosity of the aquifer, dimensionless. Krishnaswami et al. (1982) expressed R f and K in terms of adsorption and desorption rate constants, k 1 and k 2 : where i and j refer to two isotopes of the same element and Ω is the ratio of the activity of a nuclide in solution, λC, to its rate of production P: Since these authors did not consider dissolution processes, the supply rate P only includes recoil (P r ) and in situ production (λ'C') and thus, this model is only valid for radionuclides with half-lives less than 10 years (Appendix 1, Krishnaswami et al. 1982).

Precipitation and dissolution
The role of precipitation and dissolution processes becomes more apparent for longerlived radionuclides. Rate constants for precipitation and dissolution may be calculated by the mass balances of radium isotopes. Hammond et al. (1988) defined radium input to brine by dissolution processes, P w (atoms s −1 kg −1 ), by the following equation: where A is the radium activity in rocks (atoms kg −1 s −1 ), k w is the first-order rate constant for dissolution (s −1 ), and λ is the decay constant of the respective radium isotope (s −1 ).
The converse process, radium co-precipitation with minerals, S p (atoms s −1 kg −1 ), can be expressed by where k p is the first-order rate constant for precipitation (s −1 ) and C the concentration in brine (atoms kg −1 ).

Model performance
The introduced mass balance approach was modeled using Mathcad ® , a numerical software with computer algebra system (CAS) capabilities. In order to check the CAS approach, MIN3P, a multicomponent reactive transport code, was employed (Mayer et al. 2002).
MIN3P is a general-purpose flow and reactive transport code for variably saturated media providing a high degree of flexibility with respect to the definition of the reaction network. Advective-diffusive transport in the water phase and diffusive transport in the gas phase are included. Equilibrium reactions considered are aqueous complexation, gas partitioning between phases, oxidation-reduction, ion exchange, and surface complexation. The reaction network is designed to handle kinetically controlled intra-aqueous and dissolution-precipitation reactions. All reactions can be defined through databases of MINTEQA2 (Allison et al. 1991) and PHREEQC2 (Parkhurst and Appelo 1999). Table 9 lists the set of data that were used for the comparison. Input parameters correspond to the physical parameters of the Bruchsal site which are required for reactive transport modelling (Joachim et al. 1987).

Recoil mechanism
Since MIN3P treats the recoil mechanism as an intra-aqueous reaction, recoil supply rates were specified in the respective database file. Assuming that short-lived nuclides are mainly controlled by alpha-recoil, 224 Ra was used as an example to test the accuracy of the recoil term. Figure 7 illustrates the results of the comparison between MIN3P and Mathcad ® .

Kinetic approach for solid solutions
Co-precipitation of radium with barite is an important process affecting radionuclide reactive transport in rock formations. It is generally described using a solid solution model (Parkhurst and Appelo 2013). Commonly, geochemical equilibria are quantified by the law of mass action. A suitable example is given by the reaction of two components A and B with their stoichiometric constants a and b. Considering A and B as the aqueous components and AB as the solid-phase components (aA + bB ⇋ AB) leads to the equilibrium constant where brackets {} represent the activity of the components. For the condition of a homogeneous solid, its activity is assumed unity, hence Eq. (27) is simplified to K eq = {A} a {B} b , generally known as the solubility product of AB. Fig. 7 Comparison of MIN3P and CAS model for one-dimensional transport of 224 Ra along a flow line, calculated at different half-lives (t 1 = 0.5 t 1/2 ; t 2 = 2 t 1/2 ; t 3 = 45 t 1/2 ) For a solid solution, i.e., a mixture of several constituents, this simplification does not hold. The specific solid-phase activity becomes dependent on its mole fraction X i , yielding a set of concurrent equations. For simplicity, we neglect a potential non-ideality within the solid mixture, which would require the introduction of non-unity activity coefficients. It might as well be noted that often K eq is defined in a reciprocal way. In the given background the two reactions are with the associated solubility products The equilibrated solution to this set of equations is a bit more cumbersome to achieve and the reader may be referred to, e.g., Rodriguez-Galan and Prieto (2018). However, if the equilibrium can be approached in a kinetic simulation, the forward and backward reactions can be separated and make use of Lasaga's principle of detailed balancing (Lasaga 1998) to obtain a considerably more straightforward procedure: and with the kinetic rate expressions in terms of the total reaction (28a) Resulting in the forward activity product time an affinity term (1-IAP/K eq ), of which IAP is the complete ion activity product at current condition X BaSO4 /({Ba 2+ }{SO 4 2− }), and K eq is its counterpart for equilibrium conditions. Assuming that radium is always present in concentrations orders of magnitude lower compared to barium, so that the solid fraction only remains relevant for Ra 2+ , it follows That allows treating barite dissolution/precipitation as a kinetic reversible process as given in the database, independent of the radium co-precipitation process, leaving the two kinetic rate expressions (k1 + and k1−) for radium left to solve for separately, with lumped rate constants k + /k − : MIN3P's solid solution term was checked based on literature data (Grandia et al. 2008) resulting in (Ra,Ba)SO 4 equilibrium lines which have a similar slope for varying natural Ra-Ba trends (Fig. 8). The equilibrium line for the Bruchsal site shows the same trend as model-derived equilibrium lines from Grandia et al. (2008). Start and end point of the solid line is defined by the fluid-specific barium concentration and 226 Ra activity in brine. Fig. 8 Re-calculation of radium co-precipitation with barite from literature data (Grandia et al. 2008). Calculations result in (Ra,Ba)SO 4 equilibrium lines having a similar slope for varying natural Ra-Ba trends

Discussion
Modelling the physico-chemical mechanisms for radionuclide supply and removal at the geothermal site in Bruchsal are based on the following assumptions: an isotropic system is assumed in which a large conductive fracture of width w f and height h f are intersected by microfractures. These microfractures are part of the rock matrix with a very low hydraulic conductivity, and thus, water flow is assumed to take place only in the fracture.  Table 9 Site-specific input parameter used for the model comparison between Mathcad ® and MIN3P Data from Joachim et al. (1987) Parameter

Symbol Unit Value
Porosity φ -0.05 Hydraulic conductivity k f m s −1 4.3 × 10 −6 Darcy velocity v f m s −1 3.2 × 10 −8 Longitudinal dispersion α L m 0.1 The flow rate is parallel to the fracture orientation with a fracture length L f that coincides with the principal direction of groundwater flow. The fracture dimensions are assumed to be large relative to their aperture. Water-rock interactions take place in both regions, in the matrix blocks and in the fracture. However, the latter has a relatively low water-rock interaction rate due to their differences in surface area-to-fluid ratio. A schematic overview of the Bruchsal reservoir is shown in Fig. 9.
Since water-rock interaction rates occur at different timescales, their impact on U-Th series radionuclides varies depending on their half-lives. While mineral dissolution (leaching) mainly affects long-lived radionuclides, the physical process of alpha-recoil is mainly associated with short-lived nuclides (cf. Eq (17)).
Assuming that large fractures channel the geothermal brine and alpha-recoil directly into fractures is the only process adding radium into solution, the observed 224 Ra/ 228 Ra ratios should be greater than or equal to those in the rock material. However, ratios observed in the Bruchsal brine are lower (mean 224 Ra/ 228 Ra = 0.65), indicating that alpha-recoil depends on half-lives of the respective radionuclides.

Alpha-recoil supply rates and recoil efficiencies
Modeling the radium supply from recoil (P r ), Eq. (19) is applied to the site-specific input parameter. Since the recoil efficiency is hard to determine, it might be estimated from 222 Rn activity in brine. Since 222 Rn is an inert gas, it is entirely dissolved and can therefore be measured directly. Its production is from alpha-recoil of 226 Ra that is within a recoil distance of ~ 40 nm of mineral surfaces as well as from the decay of the dissolved parent 226 Ra in brine. Therefore, mass balance for 222 Rn can be expressed as follows: The fraction of 222 Rn atoms that is produced by alpha-recoil to its total fluid activity describes its recoil efficiency. Luo et al. (2000) proposed the following equation for calculating the 222 Rn recoil efficiency: assuming that ratio of alpha-recoil supply for 222 Rn and 224 Ra equals the 238 U/ 232 Th activity ratio in rocks. According to that 222 Rn recoil efficiency (ε 222 ) is ca. 23% at the Bruchsal geothermal site, suggesting that the primary source of dissolved 222 Rn is 226 Ra decay dissolved in brine (rather than the 226 Ra decay in the solid phase). Krishnaswami et al. (1982) used 222 Rn to normalize recoil efficiency. They defined the recoil efficiency of nuclide i relative to the 222 Rn efficiency as a function of (a) its position in the decay chain and (b) on the adsorption behavior of its progenitor. An example is given for the daughter-parent couple of 224 Ra-228 Ra. Since both isotopes are members of the 232 Th decay chain, 224 Ra is closely related to 228 Ra. However, while 228 Ra is generated by a single alpha decay, 224 Ra is generated by two alpha decays (cf. Fig. 1). Thus, for 224 Ra the probability of recoiling into water is significantly larger than that for 228 Ra, resulting in a (36) A 222Rn = P r,222Rn + R f ,226Ra A 226Ra .
(37) ε 222 = P r,222Rn higher recoil efficiency for 224 Ra relative to 228 Ra. Table 10 lists values of recoil efficiencies for radium as well as the recoil supply rates calculated from Eq. (19).

Diffusional flux
Since alpha-recoil supply is largest from surface areas with the highest contact areas, the recoil mechanism is most pronounced in the hydraulically inactive rock matrix. Rama and Moore (1984) suggested that diffusion along pore spaces and microfractures is believed to supply the recoiled atoms to the larger fractures where the sampled brine resides. Here, microfractures serve as diffusion pathways. The flux of radium into the large fracture fluid can be estimated from Eq. (38) as follows (Ku et al. 1992): where Φ * is the microfracture porosity, D m is the molecular diffusivity of radium, λ is the decay constant, P r is the supply rate from recoil (atoms kg −1 s −1 ). R f is the retardation factor and A is the activity of dissolved nuclide (atoms kg −1 s −1 ) with superscript (') referring to its radioactive parent. Porosity is expected to be 0.05. The molecular diffusivity for radium at 135 °C is ca. 2.6 × 10 −5 cm 2 s −1 calculated from the Stokes-Einstein relation (dynamic viscosity η = 0.285 × 10 −3 Pa s; hydrodynamic Stokes radius for radium R Ra = 3.98 Å). The term R ′ f A ′ is negligible since thorium progenitors are very insoluble in the geothermal system. The retardation factor of radium is assumed to be 1.3 estimated from the 222 Rn/ 226 Ra brine activity ratio (cf. Table 3).
Radium supply from recoil entering the fracture fluid does not only depend on diffusional flux, but also on the fracture width w f , expressed as follows: Applying Eq. (39) to radium with an estimated fracture width w f of 10 mm (as stated in the GB2 drilling report) results in a significant fractionation of the Ra isotopes (Table 12). The diffusional flux of radium considering an effective diffusion length ( φ * D m / ) is listed in Table 11. Since the effective diffusion length depends on the decay constant λ, it varies between 0.2 cm for the shortest-lived 224 Ra and several tens of centimeter for the longest-lived 226 Ra in the geothermal brine. In consequence, the discharge flux density (atoms cm −2 s −1 ) decreases with increasing decay constants limiting the 223 Ra and 224 Ra 2F Ra w f .  (Table 11). These results are in line with the statement of Rama and Moore (1984) who pointed out that migration through microfractures may restrict the input of short-lived radionuclides because the rate of diffusion through microfractures is so slow that it reduces the effect of alpha-recoil supply. However, their study focused on the very short-lived 220 Rn, and thus, one may have some doubts if it is applicable to isotopes with longer half-lives. Hammond et al. (1988), for example, suggest that radium reaches the large fractures within a few hours and thus, diffusional flux does not limit the activity of the short-lived Ra isotopes. Their study focused on the uranium and thorium series radionuclides in brines and reservoir rocks from two deep geothermal boreholes in the Salton Sea Geothermal Field (SSGF), California. From their modeling results, they postulated that the observed fluid activity of 223 Ra and 224 Ra can be explained by alpha-recoil, while only half of the 228 Ra and even less than 1% of the 226 Ra activity can be explained by alphareoil mechanisms. In consequence, the residual proportions of the 228 Ra and 226 Ra activities are contributed by weathering and leaching processes of radium from solid phases occurring on timescales comparable to the half-lives of 228 Ra and 226 Ra. Their approach is supported by an observed deficiency of 226 Ra in the SSFG reservoir section (Zukin et al. 1987).
However, in this study the situation is exactly opposite since the previous work of the authors figured out that 226 Ra is accumulated in the solids of the Bruchsal reservoir section suggesting that 226 Ra is rather removed from brine than leached from solids (Kölbel et al. 2020). 226 Ra can be removed from brine either by adsorption or by solid solution formation or both (Langmuir and Melchior 1985). Equal activities of 226 Ra and of its (unreactive) daughter 222 Rn indicate that 226 Ra is rarely adsorbed. Referring to Eq. (15) retardation factor is not only a function of k 1 and k 2 , but is also dependent on the decay constant of the respective nuclide. Should the desorption rate constant k 2 be much greater than the decay constant of 224 Ra (λ 224Ra = 2.209 × 10 −6 s −1 ), then R f,224Ra = R f,223Ra = R f,228Ra = R f,226Ra (Luo et al. 2000) applies. Otherwise short-lived radium isotopes will undergo less retardation than the long-lived 226 Ra due to their widely ranging half-lives.

Ra removal by adsorption and precipitation
Retardation factors for the short-lived radium isotopes might be calculated from the approach suggested by Krishnaswami (Eq. (22ff )). However, deviations in the P/λC Table 12 Model-derived rate constants of adsorption (k 1 ) and desorption (k 2 ); dissolution (k w ) and precipitation (k p ) of radium and the related in situ retardation factors (R r ) due to sorption processes and (R f * ) due to sorption processes as well as co-precipitation ratio from unity might rather result from 224 Ra depletion in large fractures than from adsorption. Therefore, only some general consideration about the retardation of radium could be made. Assuming that radium behaves mostly conservative ( 226 Ra/ 222 Rn = 1.3) and that the desorption rate constant is small compared to the 224 Ra decay constant which is likely because of the high radium solubility in the geothermal brine, short-lived 224 Ra should experience less retardation than 226 Ra (R f,224Ra ≠ R f,226Ra ). Hence, considering a min/max approach (R f,min = R f,224Ra = 1.0; R f,max = R f,226Ra = 1.3) with respect to the isotope half-lives, adsorption-desorption rates constants are in order of 10 −10 s −1 for k 1 and 10 −9 s −1 for k 2 (Table 12).
Previous studies of the adsorption behavior of radium in high-temperature and high-saline natural waters indicate that scavenging of radium by sorption processes might be of minor importance. Tanner (1964) proposed that during cation exchange Ra adsorption may be reduced because of the competition between radium and other alkaline earth metals for sorption sites resulting in an enrichment of radium in saline waters due to Ra-displacement from the rock surface by other cations with higher affinity to the exchange sites.
A rough estimate of the quantity of radium adsorbed on the rock surface may be calculated employing MIN3P. For calculations, an average cation exchange capacity (CEC) value of 2.00 meq/100 g was chosen since the reservoir material in Bruchsal mainly consists of quartz-dominated sandstones (CEC quartz = 0.6 meq/100 g according to Carroll 1959). Thermodynamic data for cation exchange were taken from PHREEQC2-database (Parkhurst and Appelo 1999). With respect to the chemical composition of the Bruchsal brine (Table 1), radium competes with Na, K, Mg, Sr, Ba and Cs for cation exchange sites. The results show that the adsorbed radium species vary between 10 −17 and 10 −12 meq/100 g corresponding to radium distribution coefficients K d of 0.014 to 0.016 mL g −1 .
However, these model-derived K d values may differ from the in situ distribution coefficients. Thus, the control of radium activities by adsorption cannot yet be proven without further site-specific investigations regarding the chemical behavior of radium for different environmental conditions. Should adsorption play a minor role, 226 Ra activity should be controlled by co-precipitation. Hammond et al. (1988) pointed out that the determination of first-order precipitation rate k p (s −1 ) and dissolution rate constants k w (s −1 ), respectively, may be obtained from solving simultaneously mass balance equations for radium isotopes. Combining Eq. (16) with Eqs. (18, 19, 25 and 26) and assuming C i = 0, the mass balance equation for radium at steady state can be expressed as: Simultaneous solution of Eq. (40) was performed for 226 Ra and 228 Ra, the two radium isotopes with half-lives in the order of years. The resulting rate constants are in the range of 10 −8 s −1 for both, dissolution and precipitation, whereas k p slightly exceeds k w (Table 12).
From the calculated rate constants further information can be obtained about radium kinetics. Hammond et al. (1988) argued that the mean time for radium in solution equals the mass of radium in the solid phase divided by the flux into solution expressed as follows: where M b is the mass of brine and M r is the mass of rock. Assuming a fracture porosity of 1% (since the 226 Ra accumulation in the reservoir rock was limited to fractured zones), the mean time for radium dissolution is predicted to be approximately 500 years. The mean time of radium in solution to precipitate in minerals τ p,Ra = 1/k p is estimated at circa 1 year.
The model-derived rates for dissolution (P w ) and precipitation (S p ) are listed in Table 12 demonstrating that radium will be preferred co-precipitate with minerals than leached from the solid phases. Langmuir and Melchior (1985) found that the concentrations of dissolved radium in some deep brines in north Texas were likely to be controlled by co-precipitation in sulfate minerals due to the high concentrations of sulfate and earth-alkali ions. Barite is a typical sulfate mineral incorporating radium in solid solution as [Ba,Ra]SO 4 . Both earth-alkali ions consist of an equal ionic charge and show similar ionic radii (radium = 1.52 Å, barium = 1.35 Å according to Shannon 1976). The observation by Langmuir and Melchior (1985) is confirmed by current studies dealing with the formation of Ra-bearing barite in German geothermal sites (Heberling et al. 2017;Haas-Nüesch et al. 2018).
From petrographic studies of the Bruchsal reservoir rock, it is known that barite often occurs in the reservoir section as a result of hydrothermal activities (Kölbel et al. 2020).
Zhen-Wu et al. (2016) studied barite dissolution and precipitation rates as a function of temperature and aqueous fluid composition. Their results demonstrate that barite readily achieves equilibrium with its adjacent fluid phase over a range of ionic strengths (aqueous NaCl concentrations = 0 to 1.5 molal) and in the presence of divalent metal cations (Ca, Mg and Sr) at temperatures ranging from 25 to 90 °C. They concluded that aqueous solution-barite equilibrium is broadly achieved in nature. Although reservoir temperature and molality of the NaCl geothermal brine are slightly increased compared to the experimental conditions (T = 134.7 °C; M(NaCl) = 2.1 mol kg −1 ), rate constants for k p and k w presented by Zhen-Wu et al. (2016) are in the same order of magnitude as those derived by the mass balance approach. Consequently, Ra removal from brine by co-precipitation with barite might be a possible explanation for the 226 Ra anomaly observed in the Bruchsal geothermal reservoir.

Ra supply by groundwater flow
Assuming a water recharge in the Black Forest at the Eastern main border of the URG as suggested by several numerical models of coupled heat-and fluid-flow (e.g., Clauser and Villinger 1990), meteoric water infiltrates into the Permo-Triassic aquifer of the Rhinegraben. Since this infiltrating water is low in natural occurring radionuclides, the initial concentration of radium can be assumed to be negligible (C i = 0 atoms kg −1 ). The water transit time τ b is estimated at 5000 years (assuming a flow distance of x = 5 km and a fluid velocity v f = 1 m year −1 ). Calculation of Q by applying Eq. (18) results in negative values which indicate radium loss rather than radium supply due to mass transport in groundwater. However, Q is only notable for 226 Ra (cf . Table 13), since the rate of radium loss for 228 Ra, 223 Ra and 224 Ra is ≤ 2.6 × 10 −2 atoms kg −1 s −1 .

Steady-state Ra fluid activities
At steady state, the removal rates (activity of dissolved and adsorbed radium and precipitation) are equal to inputs (radium recoil rate and production from dissolved and absorbed Th progenitors), so that (Porcelli et al. 2014): Equation (42) might be used as a control for the discussed water-rock interactions since it merges the single interaction processes which should lead to the radium fluid activities measured in the Bruchsal brine. The denominator of the fraction reflects the retardation factor R f * due to precipitation as well as adsorption and desorption. Since R f * does not differ from R f for short-lived radionuclides, R f * will increase with decreasing decay constants (Table 12). Table 13 summarizes radium production and removal rates of the water-rock interaction processes discussed. From the results it becomes obvious that for 224 Ra and 223 Ra, precipitation and dissolution processes can be neglected and so removal by decay is (42) A Ra,steady = P r + P w R f + k p = P r + P w R * f . Table 13 Model-derived production and scavenging rates (atoms kg −1 s −1 ) for radium in the Bruchsal geothermal system equal to inputs from recoil for 228 Th and 227 Th, respectively, within the solids for steadystate conditions. Consequently, alpha-recoil is the most important process for the shortlived 223 Ra and 224 Ra, to enter the fluid system. On the other hand, it is interesting to note that not only 226 Ra, but also 228 Ra is affected by dissolution-precipitation processes. This is in line with the statement of Hammond et al. (1988) who postulated that these processes occur on timescales comparable to the half-lives of 228 Ra and 226 Ra. Furthermore, the authors' previous work has shown that in the Bruchsal reservoir fractured and hydrothermally altered horizons are associated with a preferential accumulation of 226 Ra in the solid phase caused by water-rock interaction between the hot geothermal fluid and the associated solid (Kölbel et al. 2020). Studying water-rock interaction processes may therefore support the detection of productive geothermal reservoir horizons, which is one of the major challenges in geothermal exploration.

Water-rock interaction processes
Model-derived Ra steady-state activities are in good agreement with the observed Ra fluid activities (cf. Table 5) which support the applicability of the diffusional flux model. The diffusion of radium through microfractures do not only restrict the short-lived radionuclides as it is postulated by Rama and Moore (1984), but it has also a reinforcing effect on the 226 Ra activity due to its relatively high diffusional flux.

Conclusions
Even though naturally occurring radionuclides in geothermal brines can represent a challenge regarding the operational safety of a geothermal plant (there are currently some on-site research dealing with the use of inhibitors), they also can provide a supplementary tool for the characterization of geothermal reservoirs.
In this study, we investigated the behavior of U-Th series radionuclides in the brine of the Bruchsal geothermal site, located at the eastern main boundary of the Upper Rhine Graben (Germany). Permo-Triassic sedimentary rocks, affected by large-scale normal faults, host the geothermal reservoir.
Differences in chemical and physical properties result in radioactive disequilibria. Modeling the disequilibria based on radium ( 226 Ra, 228 Ra, 224 Ra, 223 Ra) enabled us to estimate rate constants of water-rock interactions. Since the daughter-parent ratio of 222 Rn/ 226 Ra is ca. 1.3, Ra retardation due to sorption processes is small resulting in adsorption-desorption rate constants in the range of 10 −10 s −1 for k 1 and 10 −9 s −1 for k 2 . Model-derived Ra distribution coefficients K d vary between 0.014 and 0.016 mL g −1 . First-order precipitation rate constant (k p = 3.4 × 10 −8 s −1 ) slightly exceeds those of dissolution (k w = 1.2 × 10 −8 s −1 ). Precipitation occurs on timescales comparable to 226 Ra