Calculating chemical speciation, pH, saturation index and calcium carbonate precipitation potential (CCPP) based on alkalinity and acidity using OpenModelica

The Calcium Carbonate Precipitation Potential (CCPP) is the amount of calcium carbonate (CaCO₃) that needs to precipitate or dissolve from water in order to reach equilibrium with CaCO₃. In equilibrium, water is saturated with CaCO₃ and the CCPP equals zero. If water is oversaturated with CaCO₃, the CCPP is positive; if it is undersaturated, the CCPP is negative (APHA/AWWA/WEF 2005).

Unlike other commonly used calcium carbonate saturation indices (such as the Saturation Index, the Relative Saturation and the Ryznar Index), the CCPP indicates both the tendency of a water to precipitate or dissolve CaCO₃ and the quantity that may be precipitated or dissolved. In drinking water treatment, knowledge about the ability of water to deposit CaCO₃ is often used in connection with the control of corrosion and scale formation in water conduits (Rossum & Merrill 1983; APHA/AWWA/WEF 2005).

In the current German Drinking Water Regulation (TrinkwV 2001), the parameter used to monitor calcium carbonate saturation is the Calcite Dissolution Capacity (Calcitlösekapazität) D_c. D_c has the same absolute value as the CCPP, but opposite sign; hence positive D_c values represent CaCO₃-undersaturated water, while negative D_c values stand for CaCO₃-oversaturated water. According to the TrinkwV, D_c must not exceed 5 mg/L at the outlet of a water works (for a description of the rationale behind this limit value see Nissing & Johannsen (2003)). Even though no minimum limit value is defined, practice shows that, preferably, water close to saturation should be distributed, as extensive precipitation of CaCO₃ may cause disruptions in the distribution network and equipment (Nissing & Johannsen 2003).

The TrinkwV further specifies that D_c is to be determined following the calculation method set out in the German standard DIN 38404-10 if the pH at the outlet is below 7.7 (water with a pH ≥7.7 is considered to meet the limit value for D_c). For this purpose, DIN 38404-10 establishes chemical principles such as the reactions and thermodynamic constants that should be at the basis of any calculation. The newest version of DIN 38404-10 (DIN 3840410:2012-12, hereafter referred to as DIN 38404-10) considers ion pairs and complex formation, including complexes with sulphate or phosphate (de Moel et al. 2013). The calculation of D_c starts out from data of the water chemistry analysis. To determine the reliability of this data, DIN 38404-10 lays down a procedure referred to as a plausibility test. To validate computer programmes, i.e. check whether their calculation of D_c is sufficiently accurate, DIN 38404-10 provides data sets of ten exemplary water samples (Stock & Johannsen 2013).

In the course of the present work, a computer programme which calculates the chemical speciation, pH, saturation index and eventually the CCPP based on the chemical principles established in DIN 38404-10 has been developed. The programme has been created in the form of an OpenModelica package. OpenModelica is an open-source modelling and simulation environment developed by the Open Source Modelica Consortium (Open Source Modelica Consortium 2017; OpenModelica can be downloaded free of charge from openmodelica.org). It is based on the object-oriented, equation-based programming language Modelica^®, devised by the Modelica Association (Modelica Association 2017). Modelica is especially useful to model and connect processes of different natural and technical domains. Moreover, in Modelica, reuse of components and evolution of models is particularly convenient (Fritzson 2015). The developed package has been conceived to become part of an OpenModelica model which combines the chemical and physical processes inside a drinking water treatment plant. Therefore, OpenModelica is considered a suitable programming environment for this project.

The programme considers all aqueous species specified in DIN 38404-10 except from Ca(OH)⁺, Mg(OH)⁺, HSO4⁻ and the ones containing phosphate. Like the standard, it applies to ‘waters, that are designated for the distribution as drinking water’ (DIN 38404-10, p. 6). Due to its intended application, the input parameters of the developed programme differ slightly from the ones in DIN 38404-10 (further details in the following section).

This paper serves readers as a guideline to develop their own programme. It follows in large part the approach presented in Eberle & Donnert (1991), but focusses on the implementation with OpenModelica. For this purpose, the six components of the developed OpenModelica package will be explained in detail, while concentrating on how the main steps of calculating the CCPP have been realised and how the different components work together. Where necessary, underlying chemical principles are elucidated. The thoroughly commented OpenModelica code of the entire package (as for example water 1) can be found in the appendix (available with the online version of this paper). From there, the code can be directly copied into the OpenModelica Connection Editor and run straightaway. In this project, OpenModelica version 1.9.6 was used.

To assess the accuracy of the developed programme, calculation results of six example waters from DIN 38404-10 will be compared to the values given in the standard and to those obtained with a programme based on a simplified model that does not consider ion pairs and complex formation.

The input parameters of the developed programme are listed in Table 1. As mentioned earlier, the programme is meant to be applied within some other larger model, for example of a drinking water treatment plant. For this reason, the developed programme starts out from the parameters ionic strength, total alkalinity and CO₂-acidity and not, like DIN 38404-10, from the pH, acid neutralising capacity with reference to pH 4.3, and base or acid neutralising capacity with reference to pH 8.2. Consequently, precision of the underlying data is also assumed to have been previously assessed. This particularly entails that the plausibility test set out in DIN 38404-10 is not within the scope of this programme. The programme also calculates the ionic strength (IOScal), but only to enable comparison with the input ionic strength (IOS). For calculation of the CCPP, IOS is used. The total concentrations of sodium, potassium, chloride and nitrate are only needed to calculate IOScal.

To simulate the effect of a specific treatment process (e.g. chlorination) on CaCO₃ precipitation or other output parameters of the programme, it is necessary to first identify how the regarded treatment process influences the programme's input parameters. This is different for each treatment process and has to be determined through measurements or models. Once the impact on the input parameters is known, the programme can be run and the process's effect on chemical speciation, pH, saturation index and CCPP evaluated.

Note that Ca(OH)⁺, Mg(OH)⁺, HSO₄⁻ as well as phosphate and its complexes are not considered within the programme. When carrying out the comparison in the end of this paper with and without considering Ca(OH)⁺, Mg(OH)⁺ and HSO₄⁻, their impact on the CCPP was found to be in the range of 0.001 mg/L, thus far beneath the accuracy limit of 0.1 mg/L demanded in DIN 38404-10. Therefore, omission of Ca(OH)⁺, Mg(OH)⁺ and HSO₄⁻ does not entail any significant losses in precision, but helps to considerably decrease code complexity and increase its readability. The neglect of phosphate is justified as its concentrations inside drinking water treatment plants are usually very small, especially compared to phosphate dosages that might be added later to inhibit corrosion. However, the developed programme can be extended to include the omitted species and reactions without any difficulties.

The function Constants calculates the concentration-based equilibrium constants Kc in mmol/L at the temperature T. As input parameters, it requires the temperature T in °C and the ionic strength IOS in mmol/L. Output parameters are the array Kc containing the eleven concentration-based equilibrium constants and the activity coefficient fH of H⁺, which will be needed in the function EqConc to calculate the pH.

All considered species, their index i and their ion size parameter can be found on the right side of Table 3. For uncharged species, lg(f(i)) is zero. The explicit activity coefficient is only calculated for H⁺, as it will be needed in the function EqConc.

The equilibrium constants K considered in the programme are listed in Table 2. They are defined based on the activities of the species involved in the respective reaction. However, for subsequent calculations it is convenient to use concentration-based constants Kc. The conversion procedure will be demonstrated exemplary for [1]. In the following, brackets [] signalize concentrations (besides identifying indices within arrays), while braces {} symbolize activities.

The function EqConc calculates the equilibrium concentration of each species in mmol/L. Required input parameters are the pH, alkalinity and acidity ⁠, the total concentrations of calcium, magnesium, sulphate, sodium, potassium, chloride and nitrate (⁠⁠, ⁠, ⁠, ⁠, ⁠, and respectively) as well as the outputs from the function Constants (⁠ and ⁠). As the input pH of EqConc is needed but unknown, it has to be determined within an iteration in the following function (function pH_Alk_Aci). Output parameters of EqConc are the array containing the 22 equilibrium concentrations, the charge factor ⁠, the residual D of the charge balance and the calculated ionic strength ⁠. Table 3 lists all components of the array ⁠.

The function begins by calculating the concentrations of H⁺ and OH⁻ from the handed over pH, activity coefficient and equilibrium constant ⁠. Then the concentrations of the main constituents (constituents are the species that do not further dissolve under the examined conditions (Eberle & Donnert 1991)) CO₃^2-, Ca²⁺, Mg²⁺ and SO₄²⁻ are calculated from the given total concentrations. The calculation will be presented taking the example of the carbonate ion CO₃²⁻.

The function pH_Alk_Aci primarily determines the equilibrium pH from the input water data. A flowchart illustrating the sequence of necessary operations is given in Figure 1. As input variables, pH_Alk_Aci requires T, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, and ⁠. The output variables are the equilibrium pH ⁠, ⁠, the equilibrium concentrations of all considered species, and ⁠. For the determination of the CCPP, the actual role of the function pH_Alk_Aci is to provide the saturation index ⁠. As described later, the other output variables do not necessarily have to be considered in the main programme; however, they are useful to determine the programme's accuracy and estimate the reliability of the resulting CCPP.

The function pH_Alk_Aci calls the two previously described functions (Constants and EqConc). As stated before, the function EqConc determines the equilibrium concentrations at any given pH. However, in the function pH_Alk_Aci, the output of EqConc is actually used to determine the equilibrium pH. This is carried out iteratively by adapting the input pH of EqConc until equilibrium state is reached with sufficient precision. Before entering the pH iteration loop, the function Constants needs to be called to provide the equilibrium constants for EqConc.

For the very first iteration step, the iteration pH is assumed to be 7. Then the function EqConc is called to determine ⁠, ⁠, D and ⁠. The returned values of and D are immediately used to evaluate Equation (20). If, for the current ⁠, the auxiliary variable is above 0, the current value of is lower than ⁠. In particular, the calculated and D are too low. Because they both increase with increasing pH (and decrease with decreasing pH), the iteration pH of the next iteration step must be raised, which is realised by setting the lower boundary pH to ⁠. If, on the other hand, is beneath 0, the higher boundary pH is altered instead. This procedure is also known as the bisection method. For the (unlikely but possible) case that is exactly 0, is the sought pH and both boundary pH values can be set to ⁠. After each iteration step, the absolute difference between upper and lower boundary pH is evaluated. If it reaches the desired precision (here 10⁻⁶), the while-loop is exited and the sought equilibrium pH is set to of the last iteration step. Note that, to enhance robustness, iteration has been limited to a pH range between 0 and 14. As this programme is developed for applications in drinking water treatment, this should not cause any problems.

The function CCPP_Alk_Aci finally determines the CCPP of the examined water. For this purpose, it needs the same input variables as the function pH_Alk_Aci (⁠⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, ⁠, and ⁠). Aside from the CCPP, it also outputs ⁠, the pH of saturation with CaCO₃. Figure 2 shows the flowchart for CCPP_Alk_Aci.

There is no way to directly calculate the CCPP. Therefore, CaCO₃ precipitation has to be simulated by iteratively withdrawing virtual amounts of CaCO₃ from the water until saturation is reached. This virtual withdrawal is realised by altering given water parameters to the way they would react to precipitation of CaCO₃. For precipitation of X mmol/L of CaCO₃, alterations are as follows:

• decreases by ⁠, as both CO₃²⁻ and Ca²⁺ contribute to the ionic strength with two times their concentration (cf. Equation (23)),

• decreases by ⁠, because CO₃²⁻ contributes to it with two times its concentration (cf. Table 1),

• increases by X, because CO₃²⁻ counteracts CO₂-acidity with one time its concentration (cf. Table 1), and

• decreases by X, as 1 mol of CaCO₃ contains the same amount of calcium.

The actual implementation of CCPP_Alk_Aci is similar to pH_Alk_Aci. Iteration is again realised through a while-loop, but this time with as iteration variable. Inside the loop, the function pH_Alk_Aci is called. Here it is crucial to alter the input parameters using as described above. The returned saturation index is then evaluated in order to determine whether saturation has been reached sufficiently. If is above 0, the altered water is still oversaturated and the current is lower than the sought CCPP. Hence, is raised for the next iteration step by setting the lower boundary CCPP to ⁠. If, on the other hand, is beneath 0, the higher boundary CCPP is adapted. For the (unlikely but possible) case that ⁠, the exact value of the CCPP is found. The loop is exited if the difference between upper and lower boundary CCPP reaches or underpasses 10⁻⁶ mmol/L. Eventually, the iteration CCPP of the last iteration step is taken as the actual CCPP (in mmol/L) of the examined water.

Note that, to enhance robustness, overall limit values for the CCPP have been set. The upper limit has been set to ⁠, meaning that only as much CaCO₃ can precipitate as calcium is contained in the water. The lower limit value has been set to −200 mg/L, as higher values are rarely found in drinking water treatment; yet it remains an arbitrary number.

The main model has to provide the input parameters listed in Table 1 in the specified units. Call the function pH_Alk_Aci in the main model to obtain the calculated pH ⁠, the sum of the concentrations of carbonic species (dissolved inorganic carbon, ⁠), the equilibrium concentrations of the examined water, the calculated ionic strength and the saturation index ⁠. This is not required to calculate the CCPP, but provides a means to determine the accuracy of the programme when testing it and to assess the reliability of its results when employing it. To obtain the CCPP (in mmol/L), the function CCPP_Alk_Aci must be called. This function also outputs the pH of calcium carbonate saturation ⁠. To obtain the CCPP in mg/L, multiply with the molar mass of CaCO₃ (here 100.09 g/mol are used). A scheme that illustrates how all described components are nested within the programme is shown in Figure 3.

To assess the accuracy of the developed programme, the and the CCPP were calculated for six example data sets of DIN 38404-10, and compared to the values given in the standard. DIN 38404-10 explicitly provides these data sets for validation of computer programmes. Data of the data sets is synthetically generated and represents analysis data of ten example waters, which mainly differ in temperature and ion content. As mentioned before, the developed programme does not consider any phosphate species. Contrary to this, the values for and of phosphate-containing waters given in the standard and required for programme execution do include these species; therefore, only the waters which do not contain phosphate (waters 1, 2, 3, 4, 7, 9) are included in this study. Table 4 summarizes relevant data of these waters. While the temperature remains constant between 10 and 15°C for all regarded waters, the ionic strength ranges widely, with a minimum at 1.041 mmol/L (water 4) and a maximum at 15.167 mmol/L (water 3). With 3.50 mmol/L and 0.70 mmol/L, water 3 contains the highest total concentrations of calcium and magnesium respectively. The highest total sulphate concentration is 1.40 mmol/L and occurs in water 7. The lowest concentrations of these three constituents are contained in water 4, with 0.15 mmol/L for calcium, 0.05 mmol/L for magnesium and 0.05 mmol/L for sulphate. Regarding their CCPP, waters 2 and 7 can be characterised as slightly, and waters 1 and 3 as moderately, CaCO₃ dissolving. Water 4 is highly dissolving, whereas water 9 is the only one with a positive CaCO₃ precipitation potential.

Values and results of the three calculation methods are presented in Table 4. The comparison reveals that CCPP values calculated with the developed programme tend to be slightly too low; yet, four of the six calculated CCPP values meet the accuracy requirement of ±0.1 mg/L set out in DIN 38404-10. With neglect of ion pairs and complex formation, i.e. with the simplified programme, no CCPP value meets the requirement. Values of comply with the required accuracy of ±0.001 for five of the six waters if calculated using the developed programme; using the simplified programme, no value lies within limits of accuracy.

Results of the simplified programme are in accordance with findings described in APHA/AWWA/WEF (2005), which states that ‘[m]odels that do not consider [ion pairs] overestimate the amount of CaCO₃ that can be precipitated and underestimate the amount of CaCO₃ that can be dissolved’ (pp. 2–34). As can be seen in Table 4, this behaviour is especially pronounced for waters with high ion concentrations (waters 3 and 9). Correspondingly, calculations for waters with low ion concentrations (waters 2 and 4) give relatively good results.

Results of the developed programme deviate from DIN 38404-10 because, contrary to what is assumed in the function CCPP_Alk_Aci, the ionic strength does not precisely decrease by four times the precipitated amount of CaCO₃. When CaCO₃ precipitates, the water's pH decreases, which affects speciation. Consequently, withdrawal of CaCO₃ not only results in lower concentrations of Ca²⁺ and CO₃²⁻, but also of other, lower charged species. Therefore, the ionic strength will generally decrease less than four times the precipitated amount of CaCO₃. To determine the exact decrease, it is necessary to additionally iterate over the ionic strength when calculating the CCPP in CCPP_Alk_Aci. However, within the additional iteration loop, has to serve as input parameter of pH_Alk_Aci. This is not desirable for this programme, as it shall also be applicable in situations in which not all constituent concentrations listed in Table 1 are known, but is assessable, for example through the electrical conductivity.

This work has shown how to programme calculation of species concentrations, pH, CaCO₃ saturation and the CCPP using OpenModelica. The advantage of using OpenModelica is that the developed package can easily be integrated into the model of, for example, a drinking water treatment plant, where it can be employed, for example, to estimate whether water leaving the plant meets the required CCPP, or to assess how much CaCO₃ dissolves if neutralizing filters containing calcium carbonate are employed to adjust the pH of the water.

The described approach adheres closely to the principles established in DIN 38404-10 and, as proven by comparison to six example waters, its results for the CCPP of phosphate-free waters nearly reach the degree of accuracy stipulated by the standard – a degree that cannot be met when neglecting ion pairs and complex formation.