*Review* **Thermodynamic Modeling of Mineral Scaling in High-Temperature and High-Pressure Aqueous Environments**

**Derek M. Hall 1,2, Serguei N. Lvov 1,2,3 and Isaac K. Gamwo 1,\***


**Abstract:** Methods of predicting mineral scale formation have evolved over the years from simple empirical fittings to sophisticated computational programs. Though best practices can now solve complex multi-phase, multi-component systems, they are largely restricted to temperatures below 300 ◦C. This review examines critical gaps in existing mineral scale modeling approaches as well as strategies to overcome them. Above 300 ◦C, the most widely used model of standard thermodynamic functions for aqueous species fails when fluid densities are below 0.7 g cm<sup>−</sup>3. This failure occurs due to the model's reliance on an empirical form of the Born equation which is unable to capture the trends observed in these high temperature, low density regimes. However, new models based on molecular solvent-solute interactions offer a pathway to overcome some of the deficiencies currently limiting high-temperature and high-pressure mineral scale predictions. Examples of the most common scale prediction methods are presented, and their advantages and disadvantages are discussed.

**Keywords:** mineral scale; high temperature; high pressure; review scaling models; aqueous species; partial molar Gibbs energy; equilibrium constants; Gibbs energy minimization

#### **1. Introduction**

The safety and reliance of many engineering applications rely on the prediction of operating scenarios that lead to severe mineral scale formation events. A classic example of the immediate and dramatic effects of scale formation is an event that occurred during early production from the Miller oil field which brought production of 30,000 barrels a day to zero in 24 h [1]. Knowledge of scaling conditions inform operating decisions that directly affect lifetime and safety [2]. As such, methods of modeling mineral solubility, which is the driving force for scale formation, date back to the early 1900s and are still being improved today. When available, mineral scale models can be used to (i) avoid mixing of incompatible brines, (ii) determine if scale inhibitors are needed and (iii) predict when scale removal maintenance is required [2,3]. Some common mineral scales are sulfate-based (Ba, Sr, Ca), oxides/hydroxide-based (Fe, Mg), and carbonate-based (Ca, Mg, Fe) [2]. However, the development of models to predict the formation of all these possibilities is an ongoing effort with many challenges due to the number of possible conditions that lead to mineral scale formation.

Determining the thermodynamic favorability of mineral scale formation is the first and most important step in combating mineral scaling. Thermodynamic favorability is governed by a mineral's solubility limit (*b*sat) and whether a particular fluid is above or below that limit. As such, knowing *b*sat is paramount to effective mineral scale control. To better predict *b*sat, the modeling of how it changes with system conditions has evolved from simple polynomials [4,5] to sophisticated thermodynamic programs [6–9] containing

**Citation:** Hall, D.M.; Lvov, S.N.; Gamwo, I.K. Thermodynamic Modeling of Mineral Scaling in High-Temperature and High-Pressure Aqueous Environments. *Liquids* **2022**, *2*, 303–317. https:// doi.org/10.3390/liquids2040018

Academic Editors: Juan Ortega Saavedra and William E. Acree, Jr.

Received: 13 August 2022 Accepted: 20 September 2022 Published: 28 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

computational algorithms and thermodynamic databases [10–12]. This progression was a necessary response to the increasing complexity and diversity of how different mineral solubility limits change with the composition, temperature and pressure of fluids observed in applications which now range from low-temperature, low-pressure (LTLP) to high-temperature, high-pressure (HTHP) conditions [2,3]. Though complex, it is vital that thermodynamic models reliably capture these trends because system properties (temperature, pressure, etc.) vary dramatically within a system which can promote scaling in different regions (see Figure 1).

**Figure 1.** (**a**) A generic demonstration of how mineral solubility limits (*b*sat) can vary with temperature for three different minerals (not to scale), (**b**) an example of how temperature and pressure gradients can promote mineral scale formation (without considering the scaling kinetics), and (**c**) an image capturing the possible consequences of mineral scale formation [13]. Figure 1c is reprinted with permission from [13]. Copyright 2022 FQE Chemicals.

The solubility limits of a mineral depend strongly on the capacity of a solution to solubilize the solute species needed to form that mineral. This capacity is called the ion activity product (*IAP*), and it is used with a solubility product (*K*sp) to assess a fluids tendency for scale formation [14,15]. This tendency is quantified with what is known as the scale saturation ratio (*SR*),

$$SR = \frac{IAP}{K\_{\rm sp}} \tag{1}$$

which is a ratio of the IAP to *K*sp for a given mineral [14,15].

For barite, a common scale in oil and gas applications, the *SR* is determined as follows:

$$SR(\text{Barite}) = \frac{a\_{\text{Ba}^{2+}(\text{aq})} a\_{\text{SO4}^{2-}(\text{aq})}}{K\_{\text{sp. BaSO\_4}(\text{s})}} \tag{2}$$

where *a*Ba<sup>2</sup>+(aq) is the activity of barium ions and *a*SO4 <sup>2</sup>−(aq) is the activity of sulfate ions within the solution. The *IAP* value represents the activities of species present in the solution, which can deviate from what is thermodynamically stable, and is affected by side reactions with other chemical species also within the solution [14]. Setting aside the non-thermodynamic contributions (e.g., kinetics of scaling), the influence of possible side reactions with additional chemical species is no trivial matter and a primary driver of the steady evolution of models to predict the solubility limit of minerals that account for multi-phase, multi-component contributions.

Many publications address the modeling of mineral solubility limits; all of which can be sub-divided into three distinct strategies that have certain advantages and disadvantages. However, the option to use a particular method depends strongly on the availability of specific datasets. Over the last 40 years, the three options most frequently used are

(i) empirical data fits of experimental data from solubility experiments, (ii) thermodynamic calculation of solubility products from thermodynamic databases, and (iii) multi-phase, multi-component computational thermodynamic programs. The first approach requires experimental solubility limit data for the system of interest, whereas the last two approaches leverage databases and publications that provide standard partial molar Gibbs energy values for the species of interest to these models. Each method can be further sub-divided to differentiate the range of sub-techniques used to obtain solubility limits of scaling minerals (see Figure 2).

**Figure 2.** Model evolution tree of mineral scale solubility models from the 1980s to today. The grey boxes are examples of the different techniques.

Though empirical fits are often the easiest to apply and the most readily available for simple systems, they have been steadily phased out from use for multi-phase, multicomponent environments.

Because a number of reviews have already been conducted on mineral scale modeling strategies, we will briefly draw attention to some of their main findings to emphasize how this review differs. A few recent publications have identified important deficiencies in scale modeling methods for the petroleum and natural gas industry up to 200 ◦C. They found that the availability of experimental data [16,17] and shortcomings in activity coefficient modeling databases are both impeding reliable mineral scale modeling efforts [16,18]. A few studies have also stressed the need to incorporate multi-phase, multicomponent features into scale modeling practices given their significant impact to several important mineral scale chemistries [16,18]. Other studies have focused on the role of scale mitigation strategies through the use of inhibitors [17,19,20]. However, none of these works discuss the underlying model limitations for the standard thermodynamic functions that underpin many of these modeling methods if extended to temperatures above 300 ◦C. Likewise, the internal consistency and viability of available standard Gibbs energy and equilibrium constant data used by many scale modeling programs are rarely discussed. Here, we provide a general overview of basic scale modeling methods, discuss overlooked technical gaps still present in commonly used standard thermodynamic functions and highlight strategies to overcome them.

#### **2. Empirical Fitting Methods**

One of the earliest methods for modeling the solubility limits of mineral scale were basic empirical fits of solubility data. Prior to the early 1980s, the partial molar standard Gibbs energy data required to determine the thermodynamic driving force for mineral scale deposition reactions were not readily accessible or known for many conditions [21]. A good example of this method is calcium carbonate mineral scale in pure water and saltwater environments [5,22]. These empirical fits were used to model solubility limits for specific systems with known compositions, temperatures, and pressures. As such, they were frequently updated as more experimental data became available; in many instances, this approach remains the best way to compile data from multiple experimental studies.

Most empirical solubility models are basic polynomials with the minimum number of parameters needed to get a satisfactory fit (i.e., a high *R*<sup>2</sup> value). These models either fit an equilibrium constant for a mineral precipitation reaction (solubility product) or a solubility limit concentration itself [4,5,22–24]. The choice of solubility product over solubility limit was often based on the desired use of the final fitting. Earlier works often focused on obtaining the solubility products, since these were not known for many mineral reactions, whereas later works focused on the solubility limits since these values are often used to validate new thermodynamic models for predicting mineral precipitation conditions.

A notable example of the simple empirical method for solubility limits was demonstrated by Krumgalz [4,23,24]. His three-part publication series analyzed over a thousand publications with mineral solubility data to formulate empirical fits of the solubility limits for alkaline and alkaline earth sulfates, chlorides and bromides at elevated temperatures. The cations examined in each publication were sodium, potassium, magnesium, calcium, strontium and barium, and the anions were chlorides, bromides, and sulfates. The temperature range varied from mineral to mineral because they were limited by the availability of experimental data. Generally, several of the minerals were fit for a range of 0 to 300 ◦C. The systems analyzed were the mineral of interest in pure water for a series of different temperatures. Each of these minerals were fit to a polynomial that had up to six of the following terms:

$$b\_{\text{stat}} = a\_0 + a\_1 t + a\_2 t^2 + a\_3 t^3 + a\_4 t^4 + a\_5 t^5 + a\_6 t^6 \tag{3}$$

where *t* is the temperature is Celsius, and *a*0–*a*<sup>6</sup> are empirical constants used to fit the experimental data. As such, this empirical model provides reliable mineral solubility limits for the water-mineral system over a given temperature range. The work by Krumgalz demonstrates both the impact of temperature on the solubility of common minerals, and how they can vary dramatically between minerals for the same set of conditions (see Figure 3).

**Figure 3.** Calcium sulfate (**a**) and barium sulfate (**b**) empirical fits of mineral solubility limits in pure water as a function of temperature [4].

Of the minerals analyzed, most fits obtained had an *R*<sup>2</sup> of 0.98 or better, with only three of the six empirical constants, but the number of constants required typically increased with the temperature range. Three to six terms are required for simple systems over small temperature ranges, it becomes increasingly difficult to accommodate additional dependences such as composition and pressure changes without a model being inundated by empirical constants.

Though less common now, this empirical approach has been used to capture the impact of more than one independent variable. Examples can be found as far back as the early 1980s, when scientists were documenting the solubility product (*K*sp) of calcium carbonate minerals in different seawater salinities for a range of temperatures [5,22]. Their goal was to model both the temperature and salinity dependences of calcium carbonate. The resulting model to capture temperature and salinity trends are as follows:

$$\log\_{10} K\_{\rm sp} = a\_0 + a\_1 T + a\_2 T^{-1} + a\_3 \log\_{10}(T) + (b\_0 + b\_1 T + b\_2 T^{-1}) S^{0.5} + c\_0 S + d\_0 S^{1.5} \tag{4}$$

where *a*0–*a*3, *b*0–*b*2, *c*<sup>0</sup> and *d*<sup>0</sup> are all the empirical constants needed to model the dependence of the solubility product on salinity (*S*), (from 5 to 44 ppt), and temperature (*T*), (from 278.15 to 313.15 K). Note that *K*sp can be converted to *b*sat through the following expression:

$$b\_{\text{sat}} = K\_{\text{sp}}^{\quad 0.5} \tag{5}$$

Using Equation (5), the impact of salinity on the solubility limit of calcite is clear (see Figure 4). Increasing salinity results in an order of magnitude increase in the solubility limit for the range of temperatures observed.

**Figure 4.** Modeled calcite solubility limits in pure water (**a**) versus in 40 ppt salinity seawater (**b**) both from 20 to 60 ◦C at atmospheric pressure [5,22].

Though good agreement can be obtained between experimental data and the resulting solubility models [4,5,22–25], the added complexity required to capture such a limited range of conditions highlights the limitations of this approach. Extending the range and properties covered by an empirical fit requires an ever-growing number of experimental data and empirical constants. In short, though this approach is the easiest model to implement, it only works for systems with available solubility data and a viable fit. Due to the strong temperature and compositional dependences of solubility, larger data sets require more empirical parameters and more empirical fits. These constraints limit the applicability of this approach to simple systems and are typically of little value to the energy industry beyond model validation.

#### **3. Solubility Product Methods**

In the 1990s, the introduction of comprehensive thermodynamic databases for minerals and aqueous species allowed for direct calculation of solubility products without fitting experimental solubility data [6,26,27]. This was accomplished with knowledge of standard Gibbs energy of formation (Δf*G*0) data for species participating in mineral scale reactions. These values can be used to calculate the standard Gibbs energy of reaction (Δr*G*0), and by extension the solubility products shown in Section 2. Unlike the empirical approaches in Section 2, this method leverages shared thermodynamic values for aqueous and solid species that can be determined by experiments not directly involving solubility measurements. This both reduces the number of empirical parameters required and expands the data available for use in mineral scale predictions. The central link between these so-called thermodynamic databases and mineral solubility limits is through the relationship between *K*sp and the Δr*G*0:

$$
\ln \left( \mathcal{K}\_{\text{SP}} \right) = -\Delta\_{\text{V}} G^0 \left( RT^{-1} \right) \tag{6}
$$

where *R* is the molar gas constant (in J mol−<sup>1</sup> K<sup>−</sup>1) and *T* is the thermodynamic temperature (in K). Despite the apparent simplicity of this approach, the underlying equations used to determine the necessary Gibbs energy data can be quite involved but they are able to capture an extraordinary range of temperatures and pressures.

One prominent example of this approach is the prediction of barite solubility between 0 and 300 ◦C in pure water using just standard Gibbs energy of formation values. For these conditions, the Δr*G*<sup>0</sup> can be calculated using the following relation which requires the partial molar Gibbs energy of formation for each aqueous species as well as the molar Gibbs energy of formation of barite at the temperature of interest:

$$
\Delta\_{\rm r}G\_{T}^{0} = \Delta\_{\rm f}G\_{\rm Ba^{2+}}^{0}{}\_{\rm (aq),T} + \Delta\_{\rm f}G\_{\rm SO^{2-}}^{0}{}\_{\rm (aq),T} - \Delta\_{\rm f}G\_{\rm BaSO\_4(s),T}^{0} \tag{7}
$$

The molar Gibbs energies of formation obtained from the SUPCRT database of minerals and aqueous species [8] required by Equation (7) results in predictions that are in excellent agreement with experimental solubility data up to 300 ◦C (see Figure 5).

**Figure 5.** Comparison between experimental data (green) [4] and model predictions (red) using Equations (5)–(7). All calculations were performed assuming negligible solution non-ideality.

Though effective for simple systems, this approach requires self-consistent thermodynamic databases which necessitate considerable work to develop. Through great efforts in the 1980s to the 1990s, a few notable databases were established and enabled solubility

product calculations over a wide range of conditions for hundreds of minerals. However, due to the considerable efforts required in building these databases, many technical gaps [6,8] have gone unaddressed in the decades since.

One of the most notable thermodynamic databases for minerals was first published by the U.S. Geological Bulletin in 1978 with revisions ending in 1984 [28]. This massive database covered thermodynamic properties of 133 metal oxides and 212 other minerals [28]. Gibbs energy values were tabulated in 100 K intervals up to 1800 K provided no phase changes were expected [28]. Temperature dependences of all mineral species were captured using a polynomial to express enthalpy at a given temperature (*H*T) relative to its enthalpy at a reference temperature (*H*298) [28]:

$$H\_{\rm T} - H\_{298} = A + BT + CT^2 + DT^{-1} + ET^{1/2} + FT^3 \tag{8}$$

where *A*, *B*, *C*, *D*, *E*, and *F* are specific to each mineral empirical constants which are fitted using available data across the temperature range. A derivation of this relationship can provide values of Δf*G*<sup>0</sup> for over 300 minerals using the following equation [28]:

$$
\Delta \mathbf{f} G\_T^0 = H\_{298} + \Delta T \left[ \frac{H\_T - H\_{298}}{T} - S\_T^0 \right] \tag{9}
$$

where *S*<sup>0</sup> *<sup>T</sup>* is the standard molar entropy of the mineral at temperature *T*. Applying this approach enabled the use of a broader range of experiments beyond solubility and mineral scale tests to fit model parameters relating to the solid phase. Shortly after the introduction of this massive mineral database, the underlying equations were also published [8] for one of the most comprehensive aqueous species thermodynamic databases.

Unlike gas and mineral phases, thermodynamic models for species in aqueous phases took much longer to create due to the complexities of solvent–solute interactions. It was not until the early 1990s that scientific literature contained a critical mass of thermodynamic properties for aqueous species at elevated temperatures and pressures required to make a comprehensive database. Much of the credit is due to a few scientists who developed the robust model underpinning most modern databases, now known as the Helgeson-Kirkham-Flowers (HKF) model, and the subsequent fitting of more than a hundred species over the next 10 years [21]. The HKF model predicts the standard partial molar Gibbs energy values of aqueous species based on the Born equation for solvation with seven empirical parameters. The HKF model works over a remarkable range of pressures (1–5000 bar) and temperatures (0–1000 ◦C) but only at relatively high density *ρ* > 0.7 g cm−<sup>3</sup> [8]. This density limitation excludes a sizable range of conditions of interest to petroleum and natural gas extraction [29] and some other important applications such as high-enthalpy supercritical geothermal technology [30] and supercritical water gasification of biomass [31] (see Figure 6).

**Figure 6.** Temperature and pressure ranges (shown in red) where the HKF model has been observed to fail due to model limitations and insufficient experimental data.

Several studies [32–37] fit available data to extend the HKF model and provide reliable predictions for multiple species that were used to form the SUPCRT database. The primary equation within the HKF model is as follows:

$$\begin{split} \Delta G\_{P,T}^{0} &= \Delta\_{\mathbb{F}} G\_{Pr,Tr}^{0} - S\_{Pr,Tr}^{0} \left( T - T\_{r} \right) - c\_{1} \Big[ T \ln \left( \frac{T}{T\_{r}} \right) - T + T\_{r} \Big] + a\_{1} (P - P\_{r}) \\ &+ a\_{2} \ln \left( \frac{\overline{\mathbb{Y}} + P}{\overline{\mathbb{Y}} + P\_{r}} \right) \left( \frac{\overline{\mathbb{Y}} + P}{\overline{\mathbb{Y}} + P\_{r}} \right) \\ &+ \Big[ a\_{3} (P - P\_{r}) + a\_{4} \ln \left( \frac{\overline{\mathbb{Y}} + P}{\overline{\mathbb{Y}} + P\_{r}} \right) \Big] \Big( \frac{1}{T - \Theta} \Big) \\ &- c\_{2} \Big[ \Big( \left( \frac{1}{T - \Theta} \right) - \left( \frac{1}{T\_{r} - \Theta} \right) \Big) \left( \frac{\Theta - \Gamma}{\Theta} \right) - \frac{T}{\Theta^{2}} \ln \left( \frac{T\_{r} (T - \Theta)}{T(T\_{r} - \Theta)} \right) \Big] \\ &+ \omega \Big( \frac{1}{\varepsilon} - 1 \right) - \omega p\_{r,Tr} \frac{1}{\left( \varepsilon\_{T,Tr} \right)} - 1 \Big) + \omega p\_{r,Tr} Y\_{Pr,Tr} (T - T\_{r}) \end{split} \tag{10}$$

where *a*1–*a*4, *c*1–*c*<sup>2</sup> and *ω* are the seven empirical constants specific to each species. In Equation (10) Θ = 228 K and Ψ = 2600 bar are constant values within the model. Ω*Pr,Tr* is the conventional born coefficient, *ε* is the relative permittivity of the solvent, and *YPr,Tr* is a constant based on of the relative permittivity of the solvent at 25 ◦C and 1 bar.

To calculate a reliable solubility product using these databases, both the empirical constants for the participating species and the underlying model equations must be valid for the conditions of interest. However, in practice one or both conditions may not be satisfied. Recent works by the authors of this review have demonstrated the failures of the HKF model in predicting the solubility of silica (see Figure 7) and barite (see Figure 8) in HTHP conditions [6,38].

**Figure 7.** Observed (•) solubility of quartz in supercritical water at *t* = 500 ◦C as a function of pressure [39] compared to NETL-PSU model [6] calculations (•) and HKF results which failed due to model limitations (-) [6]. The authors generated these plots assuming negligible activity coefficient contributions.

Models of some aqueous species were extended into the critical technical gap shown in Figure 6 by leveraging advances in molecular statistical thermodynamics [6,38]. This approach decreased the number of empirical constants from seven (aka. The HKF model) to four for ionic species or five for species with a dipole moment. It also expanded the Born equation to a more correct form that treats the solvent molecules like a particle rather than a dielectric continuum, thereby accounting for the ion-dipole and dipole–dipole interactions based on the molecular statistical theory with the mean spherical approximation (MSA). The basic thermodynamic equation in this approach for calculating the molar Gibbs energy of formation of an aqueous species is as follows [6,38]:

**Figure 8.** Comparison between HKF model (red) predictions to experimental data (green), (2700% deviations at 400 ◦C due to invalid empirical parameters and model limitations) [38]. The authors generated these plots assuming negligible activity coefficient contributions.

where Δf*G*<sup>0</sup> *Tr*,*Pr* is the molar Gibbs energy of formation of an aqueous species at *T*<sup>r</sup> = 298.15 K and *P*<sup>r</sup> = 1 bar, *a*<sup>1</sup> and *c*<sup>1</sup> are empirical parameters responsible for all short range interactions and *G<sup>k</sup> <sup>T</sup>*,*<sup>P</sup>* and *<sup>S</sup><sup>k</sup> Tr*,*Pr* are the contributions of specific molecular statistical interactions that account for hard sphere interactions [40] and electrostatic interactions [41–43]. Additionally, the standard state contributions are included to *G<sup>k</sup> <sup>T</sup>*,*<sup>P</sup>* and *<sup>S</sup><sup>k</sup> Tr*,*Pr* [44]. Within these expressions empirical diameters of the solute (σi) and solvent (σw) and dipole moment of a molecule/ion pair (*p*j) are the three parameters that are fit for an ionic species with or without a dipole [6].

The hard sphere contribution, *Gj* HS, is determined as [40]:

$$\begin{split} \frac{G\_j^{HS}}{RT} &= -\ln(1-\eta) + 3D\left(\frac{\eta}{1-\eta}\right) + 3D^2\left(\frac{\eta}{(1-\eta)^2} + \frac{\eta}{(1-\eta)} + \ln(1-\eta)\right) \\ &- D^3\left(\frac{3\eta^3 - 6\eta^2 + \beta\eta}{(1-\eta)^3} + 2\ln(1-\eta)\right) \end{split} \tag{12}$$

where *R* is the molar gas constant = 8.3145 J K−<sup>1</sup> mol<sup>−</sup>1, *η* = π*N*Aρσw3/6, ρ is the molecular density, *D* = σi/σw, *N*<sup>A</sup> is the Avogadro number = 6.0221 1023 mol<sup>−</sup>1, β = 1/(*kT*) where *k* is the Boltzmann constant = 1.3806 10−<sup>23</sup> J K−1. The mean spherical approximation (MSA) for an ion-dipole interaction was determined as [42]:

$$\frac{G\_j^{ID}}{RT} = -N\_A \varepsilon^2 z\_j^2 \frac{(1 - 1/\varepsilon)}{\sigma\_j + \sigma\_w(\beta\_6/\beta\_3)}\tag{13}$$

where *zi* is the charge number of the ionic species, *e* is the elementary charge = 1.602 10−<sup>19</sup> C, *ε* is the permittivity of the pure solvent. *ε* is related to *β*<sup>6</sup> and *β*<sup>3</sup> through the well-known Wertheim equation given as:

$$\varepsilon = \frac{\beta\_{12}^4 \beta\_3^2}{\beta\_6^6} = \frac{(1 + b\_2/12)^4 (1 + b\_2/3)^2}{(1 - b\_2/6)^6} \tag{14}$$

where *b*<sup>2</sup> is a parameter of the MSA theory. The dipole–dipole electrostatic term, *Gj DD*, is calculated as [43]:

$$\frac{G\_j^{DD}}{RT} = \frac{-8N\_A p\_j^2 (\varepsilon - 1)}{2\sigma\_w^3 \left(1 - \frac{\beta\_{12}}{\beta\_3}\right) \left(\frac{\beta\_{12}}{\beta\_6}\right)^3 + 2\varepsilon \left(\sigma\_{\dot{\gamma}} + \sigma\_w \frac{\beta\_6}{\beta\_3}\right)^3 + \left(\sigma\_{\dot{\gamma}} + \sigma\_w \frac{\beta\_{12}}{\beta\_6}\right)^3} \tag{15}$$

The Gibbs energy contribution due to the change in the standard state density is given as [6]:

$$\frac{G\_j^{SS}}{RT} = -RT\ln(\rho RT/P^\*)\tag{16}$$

where *P\** = 1 bar is the pressure of the ideal gas reference state. Lastly, the Gibbs energy that determines the difference in unit mol fraction and unit molality reference states is given as [6]:

$$\frac{G\_j^{MS}}{RT} = -RT\ln\left(M\_s/b^0\right) \tag{17}$$

where *M*<sup>s</sup> is the molar mass of the solvent = 18.015 10−<sup>3</sup> kg/mol and *b*<sup>0</sup> = 1 mol/kg is the standard molality. While this approach has been more successful than the HKF model, it has only been applied to a few dozen species and not the hundreds currently covered by the SUPCRT database.

Limitations in the HKF model and the fitted parameters can result in errors in excess of 1000% [38] when *ρ* < 0.7 g cm−<sup>3</sup> or prohibit calculations entirely [8,10]. As such, this approach to mineral scale modeling is only valid if (i) the reaction of interest is known, (ii) the database models are valid for the conditions of interest and (iii) the species of interest have available empirical constants for the conditions of interest. For these reasons, the multi-phase, multi-component systems, encountered in the applications mentioned earlier, require computational software beyond the basic solubility product methods for predicting the solubility limit of scaling minerals.

#### **4. Speciation Model Methods**

The coupled phase and thermochemistry requirements to model mineral scaling are examples of thermodynamic problems best solved using computational thermodynamics. Both empirical fits and solubility product methods are limited in their ability to predict the diverse impacts of composition-based and phase-based influences on mineral scale formation. This is in large part due to the multi-phase and multi-component nature of the scaling process. As a classical example, the impact of solution pH and CO2 partial pressure on calcium carbonate scaling are well-known [9,15,45–47], and require more complex methods to predict the extent of scaling [9,48,49]. Likewise, the advent of ion-pairs at elevated temperatures and pressures [26,50–52], complicate the identification of the correct solubility product limiting mineral solubility.

By the 1990s, speciation and phase equilibria software provided a means to leverage the thermodynamic databases with the increasing availability and computation power of computers. Phase equilibria and speciation calculations provide detailed predictions of how the solution composition changes, in addition to phase changes such as scale deposition. The impact of pH and its influence on some mineral solubility is a common example (see Figure 9).

Though many software options are available, each employs one of two methods, Gibbs energy minimization (GEM) or iteratively solving a system of equilibrium constants [10,55–57]. The most notable examples of the equilibrium constant approach are the Geochemist's workbench [58–60] and PHREEQC. It leverages the thermodynamic databases of Unitherm [10,61] and others described previously developed to solve complex speciation problems. Notable examples of the GEM approach are OLI Studio:Scalechem [12,62,63], HCh [10,61,64,65], and Thermo-Calc [57].

**Figure 9.** Comparison of speciation model predictions (black line) to experimental data (green circles) [28,53,54] for CuCl(s) over a range of HCl concentrations. The authors used mean activity coefficient data obtained for concentrated HCl solutions at 25 ◦C and 1 bar.

The equilibrium constant method requires the assembly of a systems of equations connecting species interest to potentially relevant reactions. The first step is to express each of these reactions as equilibrium constants that can be calculated using standard Gibbs energy of reactions [14]. Next, mass and charge balance equations [14] must be formulated. The resulting equations can be solved to determine the resulting speciation from known inputs.

By comparison, the GEM approach formulates an equation to calculate the total Gibbs energy of a system and minimizes the function under mass balance and charge balance constraints using Lagrangian functions [55,57]. A general description of the Gibbs energy of a system is as follows:

$$G = \sum\_{p} \Delta\_{\mathbf{f}} G\_{T\_{\nu} P\_{\mathbf{k}}}^{0} \sum\_{k} n\_{k}^{p} \tag{18}$$

where *G* is the Gibbs energy of the system and *n<sup>p</sup> <sup>k</sup>* are the molar amount of substance, *k* in phase, *p* [66]. Using the expression for the Gibbs energy of the system and equations for all necessary mass balance and charge balance constraints, a Lagrangian function, *L* can be formulated to minimize the Gibbs energy of the system as follows [13]:

$$L = G - \sum\_{i=1}^{j} \lambda\_i \Phi\_i \tag{19}$$

where *λ<sup>i</sup>* are the Lagrangian factors corresponding to each constraint function, Φ*<sup>i</sup>* [13]. By solving these series of equations, these programs provide the molar amounts of species in each phase that result in the lowest Gibbs energy for the constraints provided.

Regardless of the method, these computational thermodynamic programs permit the prediction of transitions from one predominant scaling reaction to another, which was otherwise difficult to accomplish without prior knowledge. One such example is the transition observed with barite at temperatures above 300 ◦C due to the thermodynamically favored formation of an ion pair between barium and sulfate ions (see Figure 10) [38].

**Figure 10.** The effect of forming an ion pair on the solubility of barite at high temperatures [38]. The computation program predicts that mineral solubility limit is limited by the formation of ionic species at low temperatures and the ionic pair at high temperatures. The authors generated these plots assuming negligible activity coefficient contributions.

The introduction of speciation models provides a means to resolve composition dependences on mineral scaling that limited empirical and solubility product approaches. However, it is important to note that these speciation programs are still limited by the shortcomings of the thermodynamic databases used within these programs. Therefore, technical gaps in the underlying models and limited speciation databases will impact these approaches as significantly as the solubility product method.

#### **5. Conclusions**

Multiple approaches have been developed to predict the thermodynamic favorability of mineral scale formation in aqueous systems at high-temperature, high-pressure conditions. Table 1 provides a final summary of the benefits and limitations to each method covered in this review.


**Table 1.** Overview of Mineral Scale Modeling Techniques.

Early methods were restricted to empirical fits of experimental data collected from the mineral scaling system. After the emergence of the HKF model, researchers developed thermodynamic databases able to predict the standard partial molar Gibbs energy of formation values for aqueous species over a considerable range of temperatures and pressures. Access to thermodynamic databases permitted the prediction of solubility limits for hundreds of simple mineral scale systems. Now, computational thermodynamic programs can provide solubility limits for multi-phase, multi-component systems within the confines of pre-existing database limitations.

Though many works now focus on activity coefficient models and increasing the pool of available experimental data, a third equally important factor also needs to be addressed before mineral scale modeling can be successful above 300 ◦C. As of now, almost all standard thermodynamic databases are not reliable in high temperature (>300 ◦C), low density (<0.7 g cm<sup>−</sup>3) regions. This issue stems from inherent flaws in the empirical models used to build these databases. Though models of new standard thermodynamic functions have proved to be successful, they still cover a very small number of species needed by the growing number of industrial applications in this space. Therefore, future efforts are needed to recreate these self-consistent databases which have underpinned so much recent growth in aqueous-based mineral scale modeling.

**Author Contributions:** Conceptualization, I.K.G., D.M.H., S.N.L.; writing—original draft preparation, D.M.H., I.K.G.; writing—review and editing, S.N.L., I.K.G.; visualization, D.M.H., I.K.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the U.S. Department of Energy, National Energy Technology Laboratory's ongoing research under the Advanced Offshore Research Field Work Proposal, DE-FE-FE1610260.

**Conflicts of Interest:** The authors declare no conflict of interest. This project was funded by the United States Department of Energy, National Energy Technology Laboratory, in part, through ORISE. Neither the United States Government nor any agency thereof, nor any of their employees, nor the support contractor, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

#### **References**

