Barite Scale Formation and Injectivity Loss Models for Geothermal Systems

Abstract

Barite scales in geothermal installations are a highly unwanted effect of circulating deep saline fluids. They build up in the reservoir if supersaturated fluids are re-injected, leading to irreversible loss of injectivity. A model is presented for calculating the total expected barite precipitation. To determine the related injectivity decline over time, the spatial precipitation distribution in the subsurface near the injection well is assessed by modelling barite growth kinetics in a radially diverging Darcy flow domain. Flow and reservoir properties as well as fluid chemistry are chosen to represent reservoirs subject to geothermal exploration located in the North German Basin (NGB) and the Upper Rhine Graben (URG) in Germany. Fluids encountered at similar depths are hotter in the URG, while they are more saline in the NGB. The associated scaling amount normalised to flow rate is similar for both regions. The predicted injectivity decline after 10 years, on the other hand, is far greater for the NGB ($64\%$) compared to the URG ($24\%$), due to the temperature- and salinity-dependent precipitation rate. The systems in the NGB are at higher risk. Finally, a lightweight score is developed for approximating the injectivity loss using the Damköhler number, flow rate and total barite scaling potential. This formula can be easily applied to geothermal installations without running complex reactive transport simulations.

1. Introduction

The North German Basin (NGB) and the Upper Rhine Graben (URG) are two of the main regions in Germany for geothermal exploration, exhibiting promising hydrothermal resources [1]. The success of a geothermal project greatly depends on the target reservoir’ transmissivity and temperature [2]. The high fluid mineralisation in Mesozoic sandstones [3,4,5,6,7,8], however, gives rise to challenges for their long-term utilisation [9,10]. One of which the present study focused on is barite scale formation (${\mathrm{BaSO}}_4$), originating from lowering the fluid’s state temperature and pressure during the production–injection cycle. Barite is a low-soluble scaling mineral, typically found at geothermal installations located in the NGB and the URG area that handle fluids produced from reservoirs of depths greater $2000\mathrm{m}$ [9,11,12]. Scalings are undesirable because they lower tubing diameters or reduce efficiency of heat exchangers, which leads to restoration costs [3,9,12,13]. Scale formation in the host rock at the injection site (Figure 1a), however, may constitute a more serious threat, as the related pore clogging affects the reservoir’s hydraulic properties. In case of barite scaling, this leads to irreversible injectivity loss and thus reduced overall efficiency [13,14].

Due to the prograde solubility dependency of barite, reduction of temperature along the flow path results in supersaturated conditions [16,17]. Hence, the initial reservoir fluid’s state and chemistry as well as the surface system state need to be known in advance for evaluating the scaling potential of barite for a specific geothermal site. This can be done using geochemical modelling software that applies the law of mass action together with an appropriate thermodynamic database [18]. The resulting equilibrium models yield a specific potential scaling mass based on the temperature and pressure change as well as according change in solubility. While they are easy to implement, this potential scale formation amount, however, only indicates whether precipitation can be expected and if there is a respective risk. These commonly applied equilibrium models (e.g., [19,20,21]) are insufficient in predicting the related temporal impact on injectivity because they provide no data on the distribution near the injection well. This can be achieved by using a reactive transport simulator that implements respective solute transport and a kinetic rate law. The advantage is that further site specific parameters are accounted for, such as the injection flow velocity and the precipitation rate.

Barite growth is promoted when barite-supersaturated fluids come into contact with barite in the formation rock [22]. Whether scalings grow dispersed or at specific locations has a significant impact on effective macro-scale permeability and injectivity [23,24]. It is crucial to take precipitation kinetics and flow into account to make assumptions on the expected scaling distribution in the subsurface for assessing this issue in terms of long-term utilisation of geothermal systems. Precipitation kinetics of barite have been the focus of numerous experimental studies [25,26,27,28,29,30,31,32,33], demonstrating complex dependencies on nucleation, temperature, pH, ionic strength and ion ratios. Associated pore clogging and permeability loss effects have been also studied by means of core experiments [34,35] and with regards to specific pilot sites (e.g., [13,19,36,37,38]), which highlights that barite scale formation is in fact a potential risk for geothermal systems in need of quantification and prevention measures.

Three representative cases for each geothermal region are considered in the present study. For the NGB, the geothermal plant Neustadt-Glewe (NG) was chosen, which actively produces brine from a Rhaetian sandstone aquifer [6]. Barite scalings have been observed in the filters and in the heat exchanger within the surface installations. A gradual injectivity decline over the course of many years has been attributed to this issue [13]. The geothermal site Landau (LND) was taken as a representative example for the URG. In the URG, a multi-horizon approach has been shown to be feasible, lowering exploration risk due to matrix permeability [39,40]. The well section is stretched over stratigraphical units of the Bunter sandstone, the Muschelkalk and the Permian granitic basement and also exploits a hydraulically active fault zone [39]. Furthermore, two additional cases for each region were chosen based on averaged properties, representing hypothetical sites at various depths [8]. Hereafter, they are called NGBa/NGBb and URGa/URGb, respectively. The corresponding depths and temperatures are shown in Figure 1b, where it can be seen that all sites have anomalously high temperature gradients compared to the German average [15].

This paper provides a new approach for calculating potential barite scaling amounts near the injection well and assessing the resulting impact on injectivity loss for geothermal systems in the NGB and the URG. By applying equilibrium models, the total precipitation amounts are calculated for the various sites based on the initial chemical composition of the formation fluid as well as temperature and pressure change along the flow path. One-dimensional reactive transport simulations are conducted considering radially-diverging Darcy flow and precipitation kinetics to model the constant injection of supersaturated fluids into a porous aquifer. The altered reservoir’s effective permeability and thus injectivity are assessed based on the scale formation distribution. The relevant operational parameters temperature, flow velocity and reaction rate are subjected to a sensitivity analysis in order to further provide implications for a geothermal project. Finally, a score is proposed, which aims at approximating the injectivity loss using quickly accessible parameters, applicable also to planned geothermal projects.

2. Materials and Methods

2.1. Geochemistry

2.1.1. Fluid Chemistry

Formation fluid chemistry is determined by flow path history, contact with different rock types and structures and the temperature and pressure. The amount of total dissolved solids or salinity thus can vary strongly between geothermal sites. The physicochemical parameters and chemical compositions of the considered cases are shown in

Table 1. Physicochemical parameters (upper part) and chemical composition (lower part) of the considered geothermal fluids. Measured chemical compositions are taken from given literature and converted from ($\mathrm{mg}/\mathrm{L}$) to ($\mathrm{M}$) using calculated solution densities. Reservoir chemical compositions have been calculated to achieve thermodynamic equilibrium with respect to quartz, barite, anhydrite, celesite and calcite at respective reservoir conditions, using the measured values as the basis. ${\mathrm{Cl}}^{-}$ has been additionally adjusted to achieve charge balance.

Parameter	Unit	NGBa ${}^{\mathbf{a}}$		NGBb ${}^{\mathbf{a}}$		NG ${}^{\mathbf{b},\mathbf{c}}$		URGa ${}^{\mathbf{a}}$		URGb ${}^{\mathbf{a}}$		LND ${}^{\mathbf{d}}$
		Measured	Reservoir	Measured	Reservoir	Measured	Reservoir	Measured	Reservoir	Measured	Reservoir	Measured	Reservoir
T	${}^{\circ }\mathrm{C}$	25	95	25	110	25	98	25	120	25	155	25	160
P${}^{\mathrm{e}}$	$\mathrm{MPa}$	0.1	20	0.1	30	0.1	23	0.1	20	0.1	30	0.1	30
${\rho }_{\mathrm{s}}$${}^{\mathrm{f}}$	$\frac{\mathrm{kg}}{{\mathrm{m}}^3}$	1130	1110	1190	1150	1150	1110	1050	1000	1080	1010	1070	996
$\mathrm{pH}$	−	5.60	5.33	6.00	4.91	5.20	5.64	5.80	5.66	5.50	5.28	5.15	5.41
$\mathrm{IS}$	$\mathrm{M}$	4.05	4.00	6.21	6.11	4.50	4.48	1.37	1.32	2.29	2.31	2.05	2.01
${\mathrm{Ba}}^{2+}/{\mathrm{SO}}_4^{2-}$	−	NA	3.42 $\times {10}^{-3}$	NA	2.89 $\times {10}^{-2}$	8.33 $\times {10}^{-3}$	5.47 $\times {10}^{-3}$	NA	3.19 $\times {10}^{-3}$	NA	2.64 $\times {10}^{-2}$	6.38 $\times {10}^{-2}$	3.62 $\times {10}^{-2}$
${\mathrm{K}}^{+}$	$\mathrm{M}$	2.06 $\times {10}^{-2}$	2.06 $\times {10}^{-2}$	2.87 $\times {10}^{-2}$	2.87 $\times {10}^{-2}$	2.24 $\times {10}^{-2}$	2.24 $\times {10}^{-2}$	7.88 $\times {10}^{-2}$	7.88 $\times {10}^{-2}$	1.33 $\times {10}^{-1}$	1.33 $\times {10}^{-1}$	1.06 $\times {10}^{-1}$	1.06 $\times {10}^{-1}$
${\mathrm{Na}}^{+}$		3.08	3.08	4.87	4.87	3.54	3.54	8.94 $\times {10}^{-1}$	8.94 $\times {10}^{-1}$	1.59	1.59	1.27	1.27
${\mathrm{Ca}}^{2+}$		2.11 $\times {10}^{-1}$	2.15 $\times {10}^{-1}$	3.49 $\times {10}^{-1}$	3.47 $\times {10}^{-1}$	2.33 $\times {10}^{-1}$	2.37 $\times {10}^{-1}$	1.03 $\times {10}^{-1}$	1.04 $\times {10}^{-1}$	1.82 $\times {10}^{-1}$	1.77 $\times {10}^{-1}$	1.99 $\times {10}^{-1}$	1.96 $\times {10}^{-1}$
${\mathrm{Mg}}^{2+}$		7.53 $\times {10}^{-2}$	7.53 $\times {10}^{-2}$	4.61 $\times {10}^{-2}$	4.61 $\times {10}^{-2}$	6.43 $\times {10}^{-2}$	6.43 $\times {10}^{-2}$	4.23 $\times {10}^{-3}$	4.23 $\times {10}^{-3}$	4.29 $\times {10}^{-3}$	4.29 $\times {10}^{-3}$	3.25 $\times {10}^{-3}$	3.25 $\times {10}^{-3}$
${\mathrm{Sr}}^{2+}$		5.53 $\times {10}^{-3}$	3.38 $\times {10}^{-3}$	8.57 $\times {10}^{-3}$	5.70 $\times {10}^{-3}$	5.61 $\times {10}^{-3}$	3.80 $\times {10}^{-3}$	2.58 $\times {10}^{-3}$	2.70 $\times {10}^{-3}$	4.16 $\times {10}^{-3}$	8.83 $\times {10}^{-3}$	5.09 $\times {10}^{-3}$	1.06 $\times {10}^{-2}$
${\mathrm{Ba}}^{2+}$		NA	2.71 $\times {10}^{-5}$	NA	1.30 $\times {10}^{-4}$	4.43 $\times {10}^{-5}$	3.77 $\times {10}^{-5}$	NA	1.89 $\times {10}^{-5}$	NA	9.53 $\times {10}^{-5}$	8.69 $\times {10}^{-5}$	1.10 $\times {10}^{-4}$
${\mathrm{Fe}}^{2+/3+}$		1.73 $\times {10}^{-3}$	1.73 $\times {10}^{-3}$	2.61 $\times {10}^{-3}$	2.61 $\times {10}^{-3}$	1.26 $\times {10}^{-3}$	1.26 $\times {10}^{-3}$	1.84 $\times {10}^{-3}$	1.84 $\times {10}^{-3}$	3.73 $\times {10}^{-3}$	3.73 $\times {10}^{-3}$	4.03 $\times {10}^{-4}$	4.03 $\times {10}^{-4}$
${\mathrm{Cl}}^{-}$		3.79	3.67	5.85	5.69	4.19	4.15	1.30	1.18	2.06	2.10	1.88	1.79
${\mathrm{Br}}^{-}$		4.04 $\times {10}^{-3}$	4.04 $\times {10}^{-3}$	6.73 $\times {10}^{-3}$	6.73 $\times {10}^{-3}$	5.30 $\times {10}^{-3}$	5.30 $\times {10}^{-3}$	1.29 $\times {10}^{-3}$	1.29 $\times {10}^{-3}$	2.61 $\times {10}^{-3}$	2.61 $\times {10}^{-3}$	2.84 $\times {10}^{-3}$	2.84 $\times {10}^{-3}$
${\mathrm{SO}}_4^{2-}$		5.04 $\times {10}^{-3}$	7.91 $\times {10}^{-3}$	7.58 $\times {10}^{-3}$	4.49 $\times {10}^{-3}$	5.31 $\times {10}^{-3}$	6.89 $\times {10}^{-3}$	3.74 $\times {10}^{-3}$	5.93 $\times {10}^{-3}$	3.26 $\times {10}^{-3}$	3.61 $\times {10}^{-3}$	1.36 $\times {10}^{-3}$	3.02 $\times {10}^{-3}$
${\mathrm{HCO}}_3^{-}$		3.97 $\times {10}^{-3}$	3.61 $\times {10}^{-3}$	5.97 $\times {10}^{-3}$	4.10 $\times {10}^{-3}$	7.11 $\times {10}^{-4}$	8.32 $\times {10}^{-4}$	4.21 $\times {10}^{-3}$	4.04 $\times {10}^{-3}$	7.69 $\times {10}^{-3}$	7.42 $\times {10}^{-3}$	4.01 $\times {10}^{-3}$	4.14 $\times {10}^{-3}$
${\mathrm{SiO}}_2$		NA	4.15 $\times {10}^{-4}$	NA	4.25 $\times {10}^{-4}$	NA	4.16 $\times {10}^{-4}$	NA	1.09 $\times {10}^{-3}$	NA	1.89 $\times {10}^{-3}$	2.75 $\times {10}^{-3}$	2.12 $\times {10}^{-3}$

T, temperature; P, pressure; ${\rho }_{\mathrm{s}}$, solution density; $IS$, ionic strength; NA, values are not available from sources. The sample names are North German Basin (NGB) and Upper Rhine Graben (URG) at $2000\mathrm{m}$ (a) and $3000\mathrm{m}$ (b), respectively, as well as the geothermal sites Neustadt-Glewe (NG) and Landau (LND). ${}^{\mathrm{a}}$ Wolfgramm and Seibt [8]. ${}^{\mathrm{b}}$ Naumann [6]. ${}^{\mathrm{c}}$ Kühn et al. [41]. ${}^{\mathrm{d}}$ Sanjuan et al. [7]. ${}^{\mathrm{e}}$ Down-hole pressure derived from depth (Figure 1). ${}^{\mathrm{f}}$ Calculated iteratively with PHREEQC.

. The URG has generally higher temperature gradients than the NGB (Figure 1b). Measured pH levels are all slightly acidic in the range 5.1 to 6.0.

The shown fluids are all $\mathrm{Na}-\mathrm{Cl}$ or $\mathrm{Na}-\mathrm{Ca}-\mathrm{Cl}$ dominated brines. NGB fluids are generally more saline than those from the URG. Ionic strengths for the NGB range about $4.0$–$6.2$$\mathrm{M}$ and for the URG about $1.3$–$2.4$$\mathrm{M}$). The calculated charge balance errors lie between $-2\%$ and $0.03\%$, indicating that the chemical analyses show sufficient quality. $\mathrm{Ba}$-content is seldom reported due to its low solubility.

2.1.2. Equilibrium Models and Barite Scaling Potential

The well established software phreeqc (U.S. Geological Survey) [42], version 3.6, was used to do all geochemical batch calculations. The pre-shipped database pitzer.dat was applied, which provides thermodynamic data for calculating temperature- and pressure-dependent solubility products of many common minerals, as well as coefficients for the Pitzer ion-interaction activity model [43]. It has been shown that using this database yields good results for predicting barite solubility at high ionic strengths up to $6\mathrm{M}$ [18].

The equilibrium reaction of barite in an aqueous solution is defined as:

$${\mathrm{BaSO}}_{4\left(\mathrm{solid}\right)}\rightleftharpoons {\mathrm{Ba}}_{\left(\mathrm{aq}\right)}^{2+}+{\mathrm{SO}}_{4\left(\mathrm{aq}\right)}^{2-}$$

(1)

Based on the law of mass action, the saturation state of barite in a solution is calculated with

$${\mathrm{SR}}_{\mathrm{barite}}=\frac{{\mathrm{Ba}}^{2+}·{\mathrm{SO}}_4^{2-}}{K_{\mathrm{sp},\mathrm{barite}}}·{\gamma }_{{\mathrm{Ba}}^{2+}}·{\gamma }_{S{O}_4^{2-}}$$

(2)

where ${\mathrm{SR}}_{\mathrm{barite}}$ is the saturation ratio, $K_{\mathrm{sp},\mathrm{barite}}$ is the temperature- and pressure-dependent solubility constant and ${\gamma }_i$ is the activity coefficient of the respective species. A solution is undersaturated if ${\mathrm{SR}}_{\mathrm{barite}}<1$ and supersaturated if ${\mathrm{SR}}_{\mathrm{barite}}>1$. The relationship of solubility with regards to temperature, pressure and $\mathrm{NaCl}$-equivalent ionic strength is shown in Figure 2.

Furthermore, the initial ratio of the aqueous ${\mathrm{SO}}_4^{2-}$ and ${\mathrm{Ba}}^{2+}$ concentrations constrains the total precipitation amount, which follows from Equations (1) and (2). To illustrate this, two unitless example cases are considered: (A) $a_1=1$ and $a_2=1$; and (B) $b_1=0.5$ and $b_2=2$. Both have the same initial product: $\zeta =a_1·a_2=b_1·b_2=1$, but their ratios differ. If, for both cases, the two reactants are reduced uniformly due to precipitation by 0.1, the resulting products are ${\zeta }_{\mathrm{A}}=0.9·0.9=0.81$ and ${\zeta }_{\mathrm{B}}=0.4·1.9=0.76$. The product for the case B reduces more strongly and equilibrium is reached earlier. Therefore, the stronger the initial ratio deviates from unity, the less precipitation can be expected at thermodynamic equilibrium. The resulting ion ratios of the cases are shown in Table 1.

The target Mesozoic sandstones predominantly consist of quartz. Accessory barite content is reported to exist in the host rock in both regions [44,45], with concentrations ranging from 100 to 300 $\mathrm{ppm}$ (maximum $1000\mathrm{ppm}$). It appears paragenetically next to minerals, such as anhydrite, celestite and calcite, or in fracture fillings [44]. Aqueous $\mathrm{Ba}$-concentrations are only provided for the reservoir fluids of NG and LND. Based on these compositions, the fluids appear to be close to equilibrium with respect to barite at reservoir conditions. Thus, the commonly applied assumptions is used that the fluids are in equilibrium with the mineral assemblage of the formation rock [46]. As the starting point for all consecutive equilibrium and reactive transport calculations, the concentrations reported in the literature were adjusted to achieve equilibrium with all mentioned minerals (Table 1). The ${\mathrm{Cl}}^{-}$-content was adjusted to achieve charge balance.

The maximum barite scaling amount was derived from equilibrium modelling by changing the initial P and T state with respect to the production–injection cycle (Figure 1). Re-injection temperatures usually lie in the range 45–65 ${}^{\circ }\mathrm{C}$. The system pressure usually is in the order of $1\mathrm{MPa}$ [47,48]. P increases again at the injection well, which is in the order of the reservoir’s down-hole pressure. The concentration difference between the initial and altered state represents the maximum amount of barite that can precipitate from the thermodynamic point of view.

2.1.3. Crystal Growth Kinetics

Crystal growth kinetics were considered in order to yield information on time and therefore location of scale formation. Pristine barite content in the rock matrix constitutes growth sites. Nucleation from solution can lead to an increase in the amount of active growth sites over time. This process was disregarded here, however, because the respective supersaturation ratios are assumed to be too low in the investigated cases for it to be growth determining. The reactive surface area of the available barite $S_{\mathrm{barite}}\left({\mathrm{m}}^2{\mathrm{m}}^{-3}\right)$ determines the magnitude of the reaction rate, which was approximated using

$$S_{\mathrm{barite}}={\varphi }_{\mathrm{barite}}·S·\mathrm{SF}$$

(3)

where ${\varphi }_{\mathrm{barite}}\left({\mathrm{m}}^3{\mathrm{m}}^{-3}\right)$ is the initial volume fraction of barite in the rock, $S\left({\mathrm{m}}^2{\mathrm{m}}^{-3}\right)$ is the specific inner surface area of the rock and $\mathrm{SF}$ is a dimensionless scaling factor for converting the specific into the effective reactive surface area. S was assumed to be $3·{10}^4{\mathrm{m}}^2{\mathrm{m}}^{-3}$ for the considered sandstones [2,49]. Barite was found to exist in aggregations or fracture veins [44], hence it was assumed that fluid accessibility is reduced. $\mathrm{SF}$ was therefore set to $0.05$ in order to account for the reduced volume fraction that is accessible to the pore network [50]. From measured weight fractions of the mineral assemblage, the volume fraction of barite is

$${\varphi }_{\mathrm{barite}}=(1-\varphi )\frac{\frac{w_{\mathrm{barite}}}{{\rho }_{\mathrm{barite}}}}{\frac{w_{\mathrm{barite}}}{{\rho }_{\mathrm{barite}}}+\frac{1-w_{\mathrm{barite}}}{{\rho }_{\mathrm{quartz}}}}$$

(4)

where $w_{\mathrm{barite}}\left(\mathrm{kg}{\mathrm{kg}}^{-1}\right)$ is the weight fraction of barite in the rock and ${\rho }_i\left(\mathrm{kg}{\mathrm{m}}^{-3}\right)$ is the respective density of barite ($=4480\mathrm{kg}{\mathrm{m}}^{-3}$) and quartz ($=2650\mathrm{kg}{\mathrm{m}}^{-3}$). Note that, for the sake of simplicity, only quartz and barite were assumed to be in the rock for these calculation. The weight fraction of barite in the Bunter sandstones is reported to be in the range $0.01$–$0.10$ wt% [44] in the NGB or in traces [45] in the URG, hence $0.10\mathrm{wt}\%$ was used to have a conservative estimate.

A general formulation of the reaction mechanism was applied [22]:

$$R=-k_{\mathrm{p}}(1-\mathrm{SR})$$

(5)

where $R\left(\mathrm{mol}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}\right)$ is the surface area normalised precipitation rate and $k_{\mathrm{p}}$ is the temperature and ionic strength dependent rate constant. It has been shown that it is independent from pH in the range 3–9 [30]. $k_{\mathrm{p}}$ was fitted from experimental data [32] using a first-order linear regression by minimising the averaged residuals (see Appendix A.1). The resulting empirical expression is:

$${log}_{10}{k}_{\mathrm{p}}=-\frac{2532}{(T+273)}+0.694\sqrt{\mathrm{IS}}+0.29$$

(6)

where $T{(}^{\circ }\mathrm{C})$ is the temperature and $\mathrm{IS}\left(\mathrm{M}\right)$ is the ionic strength. An extrapolation outside of these ranges can be seen in Figure 3.

2.2. Flow

2.2.1. Reservoir Hydraulics

The injectivity index $J\left({\mathrm{m}}^3{\mathrm{s}}^{-1}{\mathrm{Pa}}^{-1}\right)$ quantifies an injection well’s efficiency in terms of flow rate $Q\left({\mathrm{m}}^3{\mathrm{s}}^{-1}\right)$ and corresponding pressure build-up $\mathrm{d}P\left(\mathrm{Pa}\right)$. The applied pressure difference follows from the planned flow rate and the transmissivity of the target reservoir. It can be approximated with the Dupuit–Thiem well equation, if Darcy flow is assumed [51]:

$$Q=2\pi T\frac{s}{ln\frac{r_{\mathrm{e}}}{r_{\mathrm{w}}}}$$

(7)

where $T\left({\mathrm{m}}^2{\mathrm{s}}^{-1}\right)$ is the transmissivity, $s\left(\mathrm{m}\right)$ is the water column corresponding to the pressure build-up, $r_{\mathrm{e}}\left(\mathrm{m}\right)$ is the reach of the pressure difference and $r_{\mathrm{w}}\left(\mathrm{m}\right)$ is the well radius. Transmissivity and reach are both determined by permeability $K\left({\mathrm{m}}^2\right)$, hence solving this equation with regards to s must be done iteratively using an adequate relationship for $r_{\mathrm{e}}$ and K. However, in the following, it suffices to consider only the relative injectivity loss resulting from pore clogging. It can be shown that it approximately equals permeability loss if the reach is large compared to the well bore radius ($r_{\mathrm{e}}\gg r_{\mathrm{w}}$).

Franz et al. [52] reported that circulation rates need to be around 100 ${\mathrm{m}}^3{\mathrm{h}}^{-1}$ in order to achieve the necessary thermal output for a profitable and long-term operation. Higher reservoir temperatures compensate lower flow rates and vice versa. Further, J should be at least $50{\mathrm{m}}^3{\mathrm{h}}^{-1}{\mathrm{MPa}}^{-1}$ in order to achieve these circulation rates with realistic pressure differences. The minimum hydraulic parameters for a porous hydrothermal reservoir are given it

Table 2. Hydraulic parameters of a potential hydrothermal reservoir taken from Franz et al. [52].

${\mathit{r}}_{\mathbf{w}}\left(\mathbf{m}\right)$	$\mathit{Q}\left({\mathbf{m}}^3{\mathbf{h}}^{-1}\right)$	$\mathit{K}\left(\mathbf{mD}\right)$	$\mathit{M}\left(\mathbf{m}\right)$	$\mathit{\varphi }(-)$
0.22	100	500	20	0.2

$r_{\mathrm{w}}$ is the well radius, Q is the flow rate, K is the permeability, M is the aquifer thickness and $\varphi$ is the porosity.

. K, M and $\varphi$ constitute mean reservoir properties and were accepted to be the starting point for all cases in order to compare them with regards to barite scale formation. All reservoirs were thus simplified in the models to isotropic and homogeneous porous aquifers.

2.2.2. Reactive Transport Modelling

Flow from an injection well into a porous aquifer can be described by radially diverging flow [53,54]. If the regional hydraulic gradient is neglected and assuming homogeneous and isotropic flow properties of the aquifer, the injection well exhibits radial symmetry. Further, for fully penetrating injection wells, flow can be assumed to be planar and horizontal. Thus, a one-dimensional reactive transport model was set up to model the injection of supersaturated fluids into the aquifer. The governing equations are given in Appendix A.2.

To take barite precipitation kinetics into account, the advection–reaction equation (Equation (A5)) is modified to:

$$\frac{\partial c_i}{\partial t}=-\frac{V}{r}\frac{\partial c_i}{\partial r}-R\frac{S_{\mathrm{barite}}}{m_{\mathrm{water}}}.$$

(8)

where $c_i\left(\mathrm{M}\right)$ is the solute concentration of a respective solute species i, $t\left(\mathrm{s}\right)$ is time, $V\left({\mathrm{m}}^2{\mathrm{s}}^{-1}\right)$ is a proxy for flow ${\displaystyle \left(=Q/2\pi M\varphi \right)}$, $r\left(\mathrm{m}\right)$ is the radial distance from the well-centre, $R\left(\mathrm{mol}{\mathrm{s}}^{-1}{\mathrm{m}}^{-2}\right)$ is the barite precipitation rate (Equation (5)) and $m_{\mathrm{water}}\left(\mathrm{kg}{\mathrm{m}}^{-3}\right)$ is the solvent amount. Equation (8) is similar to Equation (A7), which was solved numerically and validated with an analytical solution (Appendix A.2). However, because the kinetic rate depends on changing activity coefficients, there is no straightforward analytical solution for this. Therefore, it was split and the resulting advective and reactive operators were solved separately by applying the sequential non-iterative approach. The upwind differences in space numerical method was used to solve the advection term ${\displaystyle \left(V\partial c_i/r\partial r\right)}$. The reaction term ($R{S}_{\mathrm{barite}}/m_{\mathrm{water}}$) was solved with phreeqc’s batch kinetic solver, which uses an implicit Runge–Kutta algorithm [42]. A regular grid with 300 nodes and a constant time step was used, which fulfilled the Courant–Friedrich–Lewy stability criteria of the upwind scheme. Temperature and pressure were set to constant injection conditions. Solute concentrations (${\mathrm{Ba}}^{2+}$ and ${\mathrm{SO}}_4^{2-}$) at the inlet boundary were taken from equilibrated, initial reservoir data. All other species were used as constant background concentration. V and $S_{\mathrm{barite}}$ were kept constant during a simulation, as feedback from reaction on flow or reactive surface area were not taken into account. The solute concentration profiles therefore reached steady-state at one point.

The shape of such steady-state concentration profiles are determined by the relationship of advection to reaction. If advection is increased, the profiles will be flattened and vice versa. A way to describe this relationship is to use the dimensionless Damköhler number $\mathrm{Da}$ [55]:

$$\mathrm{Da}\left(r\right)=\frac{r{r}_{\mathrm{c}}}{V}\frac{R}{\varphi c_{i,\mathrm{eq}}}$$

(9)

where $r_{\mathrm{c}}$ ($\mathrm{m}$) is a characteristic length set arbitrarily to $15\mathrm{m}$. It follows that $\mathrm{Da}$ is a linearly increasing function along the flow path r, since all other parameters were assumed to be constant. The steeper is this function, the shorter is the equilibrium length scale and precipitation can be expected to happen closer to the point of origin. This is illustrated in Figure 4 for various $\mathrm{Da}$-slope values. The slope of $\mathrm{Da}\left(r\right)$ can be calculated simply with:

$$m_{\mathrm{Da}}=\frac{\mathrm{dDa}}{\mathrm{d}r}=\frac{\mathrm{Da}\left(r\right)}{r}$$

(10)

or written out:

$$m_{\mathrm{Da}}=\frac{r_{\mathrm{c}}}{V}\frac{R}{\varphi c_{i,\mathrm{eq}}}=\frac{2\pi M{r}_{\mathrm{c}}}{Q}\frac{R}{c_{i,\mathrm{eq}}}$$

(11)

Barite precipitation leads to porosity decrease, which can be expressed by volume fraction:

$$\frac{\mathrm{d}{\varphi }_{\mathrm{barite}}}{\mathrm{d}t}=V_{\mathrm{m}},\mathrm{barite}R$$

(12)

where $V_{\mathrm{m},\mathrm{barite}}$ is the molar volume of barite. From this expression, the altered porosity was obtained for each domain node. Note that only quartz and barite were taken into account. To approximate the change in permeability, the widely used Kozeny–Carman relationship was applied [56,57]:

$$\frac{K_j}{K_0}=\frac{{\varphi }_j^3{(1-{\varphi }_0)}^2}{{\varphi }_0^3{(1-{\varphi }_j)}^2}$$

(13)

where ${\varphi }_0$ and ${\varphi }_j$ are the initial porosity and porosity at time step j, respectively. The effective permeability follows from a series of blocks; it was calculated using the harmonic mean of the individual blocks’ permeability [58]. Due to the radial diverging flow field, the logarithmic harmonic mean must be used:

$$K=\frac{ln\left(\frac{r_n}{r_{\mathrm{w}}}\right)}{{\sum }_{j=1}^n\frac{ln\frac{r_j}{r_{j-1}}}{K_j}}$$

(14)

where $r_j$ is the radial distance of a node j from the center of the injection well and $r_{\mathrm{w}}$ is the well radius. $r_n$ was chosen so as to capture the saturation length scale.

The order of steps used for assessing the temporal permeability loss is as follows:

• Calculate reactive surface area (Equation (3)).
• Calculate rate constant (Equation (6)).
• Calculate flow constant V.
• Evaluate Damköhler number along the r-axis (Equation (9)) and the corresponding slope (Equation (11)).
• Solve ARE (Equation (8)) with numerical simulations until solutes reach steady-state.
• Extrapolate porosity change at steady-state for each node over ten years.
• Evaluate permeability loss using both the porosity–permeability relationship (Equation (13)) and the effective permeability expression for the radial diverging flow field (Equation (14)).

Each geothermal sample case was evaluated this way, representing the respective base scenario by using values described in the previous sections (Table 1 and Table 2). To illustrate the effects of various parameters on the simulation results, additional scenarios were examined as in a one-at-a-time sensitivity analysis (

Table 3. Varied parameters in the respective scenarios for a one-at-a-time sensitivity analysis. Decreasing Q and $w_{\mathrm{barite}}$ corresponds to decreasing flow velocity and precipitation rate, respectively.

Scenario	Parameter	Value	Unit
$T+10{}^{\circ }\mathrm{C}$	T	65	${}^{\circ }\mathrm{C}$
$T-10{}^{\circ }\mathrm{C}$	T	45	${}^{\circ }\mathrm{C}$
$Q/2$	Q	50	${\mathrm{m}}^3{\mathrm{h}}^{-1}$
$R/10$	$w_{\mathrm{barite}}$	0.01	−

).

3. Results

3.1. Equilibrium Models

The scaling potential results from equilibrium calculations are shown in Figure 5. The whole range from the respective reservoir temperature down to $25{}^{\circ }\mathrm{C}$ illustrates the effect of temperature reduction, expressed as $\mathrm{d}T=T_{\mathrm{res}}-T_{\mathrm{i}}$. Furthermore, the influence of pressure reduction $\mathrm{d}P$ from reservoir to surface condition is shown.

Figure 5a presents the respective barite saturation associated to $\mathrm{d}T$ and $\mathrm{d}P$. The curves show an exponentially increasing trend for all cases with increasing $\mathrm{d}T$. The URG cases (grey tones) generally exhibit higher saturation ratios than the NGB cases (red tones) at equal depths. The considered pressure reduction additionally increases the saturation, however not as strongly as the accompanying temperature reduction. With regards to the flow path in geothermal installations, saturation increases during the passageway through the production well and the heat exchanger due to temperature and pressure reduction, respectively. Supersaturation reaches its highest magnitude up to the point where the fluid is pressurised again and re-injected through the injection well. Assuming a temperature and pressure reduction down to $55{}^{\circ }\mathrm{C}$ and $1\mathrm{MPa}$, respectively, $\mathrm{SR}$ values for the various cases lie between $3.2$ and $7.1$. If these values overstep the supersaturation threshold, nucleation can be expected in the surface installation. $\mathrm{SR}$ values representative for the injection location are between $2.7$ and $5.3$, which are lower again due to the increased fluid pressure.

Figure 5b presents the associated total amount of barite that can precipitate from a cubic metre of produced formation fluid, which has been subjected to change in temperature and pressure. The cases representing greater depths (NGBb, URGb and LND) have a higher scaling potential than the shallower cases (NGBa, URGa and NG). With regards to fluid injection, values lie between 74 and 88 $\mathrm{mmol}{\mathrm{m}}^{-3}$ and 12 and 23 $\mathrm{mmol}{\mathrm{m}}^{-3}$ for the deeper and shallower cases, respectively. The results at injection conditions are also shown in

Table 4. Summary of results of the equilibrium calculations and the reactive transport simulations for all considered geothermal cases and scenarios. Further, the developed empirical scaling score $X_{\mathrm{score}}$ is shown (Equation (15)).

	Scenario	${\mathit{n}}_{\mathbf{barite}}$	${\mathit{m}}_{\mathbf{Da}}$	Loss	${\mathit{X}}_{\mathbf{score}}$
Case		$\left(\frac{\mathbf{mmol}}{{\mathbf{m}}^3}\right)$	$\left(\frac{1}{\mathbf{m}}\right)$	$\left(\frac{1}{\mathbf{year}}\right)$	$\left({10}^{-4}\frac{\mathbf{mol}}{\mathbf{m}\mathbf{year}}\right)$
NGBa	Base	16	4.4	0.018	6.1
NGBa	$T+10{}^{\circ }C$	16	5.6	0.018	7.9
NGBa	$T-10{}^{\circ }C$	16	3.3	0.016	4.7
NGBa	$Q/2$	16	8.7	0.016	6.1
NGBa	$R/10$	16	0.43	0.0026	0.61
NGBb	Base	87	2.8	0.064	22
NGBb	$T+10{}^{\circ }C$	87	3.6	0.069	27
NGBb	$T-10{}^{\circ }C$	87	2.1	0.056	16
NGBb	$Q/2$	87	5.6	0.06	22
NGBb	$R/10$	87	0.28	0.01	2.1
NG	Base	23	4.1	0.024	8.3
NG	$T+10{}^{\circ }C$	23	5.3	0.025	11
NG	$T-10{}^{\circ }C$	23	3.2	0.022	6.4
NG	$Q/2$	23	8.2	0.022	8.3
NG	$R/10$	23	0.41	0.0036	0.83
URGa	Base	12	1.6	0.0057	1.7
URGa	$T+10{}^{\circ }C$	12	2.3	0.0066	2.4
URGa	$T-10{}^{\circ }C$	12	1.2	0.0048	1.2
URGa	$Q/2$	12	3.3	0.0051	1.7
URGa	$R/10$	12	0.16	0.0008	0.17
URGb	Base	74	1.1	0.024	6.9
URGb	$T+10{}^{\circ }C$	74	1.4	0.029	9.2
URGb	$T-10{}^{\circ }C$	74	0.79	0.02	5.1
URGb	$Q/2$	74	2.1	0.022	6.9
URGb	$R/10$	74	0.11	0.0034	0.69
LND	Base	85	0.8	0.022	6
LND	$T+10{}^{\circ }C$	85	1	0.026	7.8
LND	$T-10{}^{\circ }C$	85	0.58	0.018	4.3
LND	$Q/2$	85	1.6	0.02	6
LND	$R/10$	85	0.08	0.003	0.6

T is the injection temperature, Q is the flow rate, R is the precipitation rate, $n_{\mathrm{barite}}$ is the total precipitation potential, $m_{\mathrm{Da}}$ is the slope of the Damköhler number along the r-axis and Loss represents the relative, effective permeability decrease per year, as in $1-K/K_0$.

. Lower injection temperatures result in more scaling potential, although the curve flattens off, meaning that it does not grow linearly. This further indicates that the potential scaling amount is not linearly correlated to the previously reported exponentially increasing saturation ratios. Generally, the higher the fluid’s ionic strength, the more scaling can be expected. The same applies to the degree of temperature and pressure reduction, although the drop from reservoir to system pressure has a generally smaller effect on precipitation than temperature. The NGB cases tend to exhibit higher values, although their $\mathrm{d}T$ values to reach injection temperature are lower. For instance, temperature reduction for the NGBb case is close to half compared to the URGb case, however scaling potential is almost $20\%$ higher for the NGBb case.

The initial ion ratios of aqueous ${\mathrm{SO}}_4^{2-}$ and ${\mathrm{Ba}}^{2+}$ are shown in Table 1. ${\mathrm{SO}}_4^{2-}$ concentrations are generally higher, between 30 and 300 times the concentration of ${\mathrm{Ba}}^{2+}$. The ratios of the deeper cases are about one order closer to unity than the shallower cases. Cases with ion ratios closer to unity are reflected in steeper curve slopes in Figure 5b.

3.2. Reactive Transport Models

The Damköhler number increases linearly for the considered flow model, originating from zero along the r-axis. The associated slopes $m_{\mathrm{Da}}$ (Equation (10)) for all scenarios are shown in Table 4. The higher is $m_{\mathrm{Da}}$, the faster is the advection time decreases along the r-axis compared to the reaction rate, and thus the more precipitation will occur more concentrated close to the injection well. All URG cases have comparably low slope values, whereas the NGBa and NG cases have the highest values. Reducing the flow rate increases while reducing the reaction rate reduces $m_{\mathrm{Da}}$ proportionally. In accordance with the kinetic rate expression (Equation (6)), increasing temperature increases $m_{\mathrm{Da}}$.

In the following, the reactive transport simulation results are presented. In Figure 6 and Figure 7, the NGB and URG cases are shown, respectively. On the left-hand side, the temporal porosity changes at steady-state are plotted against the r-axis. Porosity changes result from increase of barite volume fraction. Note that the porosity change plots are cut off at ${10}^{-4}$ per year, as lower values are assumed to be negligible with regards to their impact on permeability during the life time of a geothermal installation. This corresponds to a porosity loss of $-0.1\%$ over the course of ten years. All curves show flipped parabola shapes in the semi-log plot, with the porosity change maximum close to the origin, i.e., the outflow of the injection well. The NGBb case generally exhibits the highest porosity changes with a maximum of $0.7\%$ porosity loss per year. The URGa case has the least loss of more than one order of magnitude less than that of NGBb. It can be seen that the URG cases have flatter curves and broader widths of significant precipitation along the axis. Precipitation is concentrated closer to the well for the NGB cases. They have ranges of about $4.5$–$6.5$ $\mathrm{m}$, compared to $5.5$–$10.5$ $\mathrm{m}$ for the URG cases.

Changing the injection temperature within in the considered ranges of the sensitivity analysis only has a small effect compared to the base case. A temperature increase accelerates precipitation and vice versa. As such, lower injection temperature increases the reach and also lowers the maximum porosity decline close to the well. Smaller flow rates generally result in less precipitation in total, which can be anticipated by the total area under the curves; compared to the base case, the maximum porosity changes at the inlet still has the same order of magnitude, however the reach is reduced by about $30\%$. Lower reaction rates (R) flatten the curves quite significantly, meaning precipitation is more distributed along the flow path. If the magnitude of R is one order lower compared to the base case, the porosity change peak is also one order of magnitude lower and the reach is increased by about $50\%$.

The associated temporal permeability ratios ($K/K_0$) over the course of ten years are shown on the right-hand side of Figure 6 and Figure 7. Notably, all curves illustrate a linear decline of permeability over time. The related injectivity losses per year ($1-K/K_0$) are shown in Table 4. In the base case, the NGBb case exhibits the strongest permeability decline by over $6\%$ per year. It slightly deviates from a linear curve towards the end of the considered time. The URGa case shows the lowest permeability decline of just below $0.6\%$ per year; the others have between $1.8\%$ and $2.4\%$ loss per year. With regards to the sensitivity analysis scenarios, change in reactivity has the strongest impact on the permeability decline. If the reaction rate is reduced by one order of magnitude compared to the base case, the resulting permeability loss after 10 years is also reduced by almost one order of magnitude. For the shallower cases (NGBa, URGa and NG), the other sensitivity analysis scenarios with regards to temperature and flow rate show only negligible deviations from the base case. For the deeper cases (NGBb, URGb and LND), increasing temperature (lower $\mathrm{d}T$) accelerates permeability loss and decreasing temperature (higher $\mathrm{d}T$) slows the decline down. Decreasing the flow rate reduces the permeability loss. However, this effect is small, about $10\%$ deviation from the base case.

The final effective permeability changes after ten years for all investigated cases and scenarios are summarised in Figure 8. It can be seen that the reservoirs of greater depth generally are affected more strongly by permeability losses. The NGB cases are generally affected more strongly than the URG cases; NGBb has by far the highest permeability loss ($64\%$). Further, NG is similarly affected than URGb and LND, although it is not as deep; all exhibit a loss of approximately $24\%$. URGa has the least loss of just below $6\%$. Again, for all shallower cases (NG, NGBa and URGa), varying injection temperature or flow rate only has a little to no effect compared to the base case. There is a noticeable effect, however, with regards to injection temperature variation for the deep cases. As such, decreasing temperature reduces the permeability loss, whereas increasing temperature results in more loss. If reactivity is reduced compared to the base case, permeability loss is significantly less. Especially for the NGBb case, loss after ten years drops down from $64\%$ to only $10\%$ in this scenario. The loss of the URGa case goes down to virtually zero in this scenario.

4. Discussion

Total barite scaling potential is determined by fluid composition, temperature reduction, pressure reduction and ion ratio of ${\mathrm{Ba}}^{2+}$ and ${\mathrm{SO}}_4^{2-}$. For the presented geothermal cases in the NGB and the URG, formation fluid temperature and salinity increase with greater reservoir depth. The deeper cases also showed ion ratios closer to unity. These factors all increase the precipitation potential and thus increase the scaling risk for deeper reservoirs. Furthermore, for cases at similar depths, fluid salinity is higher in the NGB and temperature is higher in the URG, but total scaling potential is in the same order of magnitude. As the corresponding saturation ratios were significantly higher for the URG cases, this illustrates that these values are not linearly correlated. Saturation ratios can only be taken as an indication of whether precipitation can be expected or not. Scaling amounts need to be investigated in detail for a specific location, taking the mentioned factors into consideration.

An irreversible injectivity decline of about $35\%$ after 16 years of injection has been reported at the geothermal site Neustadt-Glewe [13]. The authors attributed this mainly to formation of sparingly soluble scales in the reservoir, such as barite, celestite and various sulfides [9,13]. This order of magnitude is in fact in accordance with calculated injectivity losses, even though only barite scaling was taken into consideration in the present study. In this regard, some conservative assumptions were made, to the effect that the upper range of risk associated to barite scaling was assessed. For instance, it was assumed that barite growth does not happen until the re-injected fluid comes into contact with active growth sites (solid barite) in the formation rock. Precipitation rate and injectivity loss are further increased. Formation of additional active growth sites can occur through nucleation, increasing the precipitation and injectivity loss rate further. Whether this effect is significant depends on site characteristics such as mineralogy of the in-situ rock. This process was disregarded in the models, hence the real precipitation rate will be underestimated to some degree. For the considered cases, nucleation was not considered to be a growth determining step, as the highest saturation ratios were presumed to be too low for this [59,60]. Nuclei formation can be promoted, however, by longer shut-in periods [10,11,36]. As a consequence, scale formation sets in prior to the fluid reaching the formation rock, thus not affecting the injectivity as strongly, but perhaps with other unwanted impediments [10]. Furthermore, re-injected fluids heat-up again gradually in the reservoir, depending on flow rate and heat transfer, which has not been taken into consideration. This increases barite solubility again and thus reduces scaling risk, more so if flow rate is reduced. In light of the simulated scaling reach of below $10\mathrm{m}$, however, this effect appears to be negligible. On the other hand, injection pressure must be increased to maintain injection rates if permeability decreases. This could pose problems, as this increases the chance for loose particles to redistribute and clog pores. This process, however, is hard to quantify and also was not considered.

Supersaturation and kinetic rate both depend on the fluid’s salinity and temperature. An increase in salinity therefore increases the precipitation rate two-fold. The relationship with regards to temperature is different: supersaturation increases further with temperature reduction (increased $\mathrm{d}T$), whereas the rate constant is proportional to the absolute temperature. For quartz scaling in high-enthalpy systems, something similar was shown by Pandey et al. [61]. Reducing temperature results in a counter effect of higher supersaturation, but lower kinetic rate constant. Temperature variations as part of the sensitivity analysis had no significant impact on calculated injectivity loss, at least in the considered $\pm 10{}^{\circ }\mathrm{C}$ range. This explains why the NGB cases are generally affected more strongly.

Scaling potential needs to be put into perspective with regards to distribution along the flow path in order to assess implications for system longevity. Due to the radial diverging flow, large spatial distribution of scale formation means less effective permeability loss. This is promoted by slower reaction rates and higher flow velocities. The URG cases generally showed widespread distribution along the flow path, i.e., flatter precipitation curves, attributable to lower reaction rates. An important point is that equal hydraulic properties were assumed for all cases, in order to be able to compare them with regards to fluid chemistry and precipitation kinetics. While the model assumptions of radial diverging and planar flow are reasonable for homogeneous and isotropic porous aquifers with fully penetrating injection wells, they are simplifying for partially penetrating wells and especially fractured aquifers. The former have spherical flow components, thus this model treatment overestimates flow velocities and underestimates the permeability loss in the near-vicinity of the injection well for these cases. Projects in the URG rely on multi-horizontal approaches in order to minimise exploration risk [39]. Therefore, there are sections with slow flow through the porous matrix, but also sections with increased permeability due to fractures. Fracture permeability is characterised by preferential flow with increased flow velocities and also decreased water–rock contact, i.e., less effective reactive surface area. Both factors hypothetically increase the scaling distribution in the formation rock, reducing the scaling risk for the URG even further compared to the NGB.

Scaling distribution patterns in the subsurface can be described by fluid flux, total scaling potential and Damköhler number. The latter relates the respective magnitudes of advection and precipitation kinetics. Assuming that rock reactivity is homogeneous on a large scale, the Damköhler number increases linearly along the r-axis due to the radially diverging flow. The steeper is the slope ($m_{\mathrm{Da}}$), the closer isthe scaling distribution to the injection well. Further, the volumetric scaling potential ($n_{\mathrm{barite}}$) as well as fluid flux (Q) are also necessary to determine scaling in the subsurface with regards to injectivity decline, as these quantify the total amount that can precipitate from solution. In essence, these three factors provide insights into distribution and intensity of scaling in the subsurface. If they are simply lumped together, the following scaling score is derived:

$$X_{\mathrm{score}}=m_{\mathrm{Da}}·n_{\mathrm{barite}}·Q$$

(15)

The resulting scores for the respective cases and scenarios are provided in Table 4. They qualitatively suggest which cases’ injectivity will be affected more than others and therefore generate a ranking fit. Although this score only yields an approximation, it can nevertheless be used as a quick comparative value, without having to run elaborate reactive transport simulations. Furthermore, this scaling score is correlated with the previously calculated injectivity losses. For instance, the NGBb case has the largest values, while URGa has the lowest, which corresponds closely to the simulation results. If plotted against each other, a clear linear correlation can be seen (Figure 9). This is a valuable insight, since the calculated injectivity losses result from multiple non-linear considerations: (I) steady-state reactive transport simulations; (II) porosity–permeability relationship; and (III) effective permeability approximation. By calibrating the score with the reactive transport simulation results, the obtained linear correlation represents a lightweight score for approximating the temporal injectivity loss associated to barite scaling:

$$\mathrm{Loss}(\%{\mathrm{year}}^{-1})=2.89·{10}^{-5}·X_{\mathrm{score}}\left(\mathrm{mol}{\mathrm{m}}^{-1}{\mathrm{year}}^{-1}\right)$$

(16)

It is easily applicable to new geothermal installations and may be calibrated further if additional data becomes available in this regard. The overall presented approach specifically treats barite as the sole scale formation agent. It can be adapted to make respective predictions for similar formation reactions of minerals exhibiting prograde solubility, for example silica or other sulfates. Though it is explicitly pointed out that the respective reaction mechanism needs detailed consideration.

5. Conclusions

Two model concepts were presented to approximate barite scaling formation in geothermal systems of the North German Basin and Upper Rhine Regions regions: an equilibrium model approach and a transport model coupled with precipitation kinetics. It was shown that temperature and pressure reduction during the production–injection cycle results in supersaturated conditions for barite in all cases, which is accountable for scaling. Equilibrium models were used to calculate the total potential scaling amount, which gives a first indication on the related risk for a long-term operation. This scaling potential increases proportionally to the imposed degree of temperature and pressure reduction dependent on the respective geothermal system management, as well as to formation fluid salinity. These parameters are generally correlated with reservoir depth. Fluids encountered at similar depths are hotter in the URG, while they are more saline in the NGB. The scaling potential is similarly high for both regions, while deeper reservoirs tend to be affected more strongly.

A comprehensive assessment of scaling risk needs to include the respective scaling location and distribution along the flow path in order to quantify the accompanied injectivity decline. From reactive transport simulations, information on both the scaling distribution in the subsurface and the related injectivity loss was obtained. Precipitation kinetics are taken into account, which also depend on temperature and salinity, similarly to the total scaling potential. Injection temperatures are usually in the same order for different geothermal installations, thus the corresponding temperature reduction ($\mathrm{d}T$) varies. The barite precipitation rate is higher for the NGB cases due to their higher fluid salinities. Thus, scaling will preferentially happen closer to the injection well and damage reservoir permeability more severely. Therefore, the NGB cases are generally at higher risk with regards to injectivity losses, while the shallow URG case showed almost no losses. A sensitivity analysis showed that varying temperature within a $10{}^{\circ }\mathrm{C}$ margin, as well as significantly reducing the flow rate had negligible effects on injectivity loss. The kinetic rate, on the other hand, exhibited a strong sensitivity.

A scaling score was developed, which takes the total scaling potential, the Damköhler number and the flow rate into account. It correlates strongly with the results of the reactive transport simulations and may be calibrated with further data. It is easily applicable in order to get an indication on the accompanied scaling risk for a specific geothermal location, without having to run elaborate reactive transport simulations. The presented approach can be adapted to make scale formation and injectivity loss predictions for mineral formation reactions similar to that of barite.