Abstract
Of the most commonly used calcium carbonate (CaCO3) saturation indices, the Calcium Carbonate Precipitation Potential (CCPP) most accurately describes CaCO3 precipitation and dissolution. This paper illustrates how to write an OpenModelica programme that calculates the CCPP and other key water parameters starting from total alkalinity and CO2-acidity. It describes the calculation steps of all components of the developed OpenModelica package and explains how the components work together. Underlying chemical principles follow the German standard DIN 38404-10 and are presented along with the OpenModelica code. The developed programme considers 14 different dissolved species, including complex species of carbonate, calcium, magnesium and sulphate. Data of the six phosphate-free example waters from DIN 38404-10 is compared to the results of the developed programme and the results of a programme which neglects complexation. The developed programme reaches required accuracy for four of the six CCPP values; with the other programme, requirements cannot be met. By using OpenModelica, the developed package can be readily integrated into a Modelica model of, for example, a drinking water treatment plant. There, it can help to monitor legal and technical requirements, or to assess the dissolution capacity if neutralizing filters containing CaCO3 are employed for pH adjustment.
INTRODUCTION
Calcium carbonate precipitation potential
The Calcium Carbonate Precipitation Potential (CCPP) is the amount of calcium carbonate (CaCO3) that needs to precipitate or dissolve from water in order to reach equilibrium with CaCO3. In equilibrium, water is saturated with CaCO3 and the CCPP equals zero. If water is oversaturated with CaCO3, 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 CaCO3 and the quantity that may be precipitated or dissolved. In drinking water treatment, knowledge about the ability of water to deposit CaCO3 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) Dc. Dc has the same absolute value as the CCPP, but opposite sign; hence positive Dc values represent CaCO3-undersaturated water, while negative Dc values stand for CaCO3-oversaturated water. According to the TrinkwV, Dc 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 CaCO3 may cause disruptions in the distribution network and equipment (Nissing & Johannsen 2003).
The TrinkwV further specifies that Dc 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 Dc). 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 Dc 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 Dc is sufficiently accurate, DIN 38404-10 provides data sets of ten exemplary water samples (Stock & Johannsen 2013).
The developed programme
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.
METHODS
Input parameters of the developed programme and differences to DIN 38404-10
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 CO2-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.
Input parameters of the developed programme
Parameter . | Description . | Unit . |
---|---|---|
![]() | Water temperature | °C |
![]() | Ionic strength at temperature ![]() | mmol/L |
![]() | Total alkalinity: ![]() | mmol/L |
![]() | CO2-acidity:![]() | mmol/L |
![]() | Total concentrations of calcium, magnesium, sulphate, sodium, potassium, chloride or nitrate, respectively | mmol/L |
Parameter . | Description . | Unit . |
---|---|---|
![]() | Water temperature | °C |
![]() | Ionic strength at temperature ![]() | mmol/L |
![]() | Total alkalinity: ![]() | mmol/L |
![]() | CO2-acidity:![]() | mmol/L |
![]() | Total concentrations of calcium, magnesium, sulphate, sodium, potassium, chloride or nitrate, respectively | mmol/L |
To simulate the effect of a specific treatment process (e.g. chlorination) on CaCO3 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)+, HSO4− 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 HSO4−, 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 HSO4− 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.
Determining concentration-based equilibrium constants (function Constants)
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.
Definition of the activity-based equilibrium constants K, their decimal logarithm at 25°C, the reaction enthalpy
in J/mol and the heat capacity
of the reaction in J/(mol·K); values are adapted from DIN 38404-10; the indices listed correspond to the indices of
in the programme's array
Index . | Activity-based equilibrium constant . | lg(K0) . | ΔH . | Cp . |
---|---|---|---|---|
1 | ![]() | −6.356 | 9,250 | −330 |
2 | ![]() | −10.329 | 14,950 | −290 |
3 | ![]() | −1.121 | −7,950 | 0 |
4 | ![]() | −3.200 | −16,000 | −430 |
5 | ![]() | −2.310 | −7,600 | −170 |
6 | ![]() | −1.068 | −7,250 | 0 |
7 | ![]() | −2.947 | −13,000 | −70 |
8 | ![]() | −2.263 | −18,000 | −10 |
9 | ![]() | −13.996 | 56,532 | −197 |
10 | ![]() | −8.481 | −10,000 | −250 |
11 | ![]() | −18.156 | 4,150 | −460 |
Index . | Activity-based equilibrium constant . | lg(K0) . | ΔH . | Cp . |
---|---|---|---|---|
1 | ![]() | −6.356 | 9,250 | −330 |
2 | ![]() | −10.329 | 14,950 | −290 |
3 | ![]() | −1.121 | −7,950 | 0 |
4 | ![]() | −3.200 | −16,000 | −430 |
5 | ![]() | −2.310 | −7,600 | −170 |
6 | ![]() | −1.068 | −7,250 | 0 |
7 | ![]() | −2.947 | −13,000 | −70 |
8 | ![]() | −2.263 | −18,000 | −10 |
9 | ![]() | −13.996 | 56,532 | −197 |
10 | ![]() | −8.481 | −10,000 | −250 |
11 | ![]() | −18.156 | 4,150 | −460 |
Components of the array ; values for
are taken from DIN 38404-10 and Eberle & Donnert (1991)
Total concentrations . | Species concentrations . | |||
---|---|---|---|---|
Index . | Component . | Index . | Component . | Ion size parameter ![]() |
1 | tC | 9 | H+ | 6.8 |
2 | tCa | 10 | OH− | 4.2 |
3 | tMg | 11 | CO2 | 3.0 |
4 | tSO4 | 12 | HCO3− | 1.5 |
5 | tNa | 13 | CO32− | 4.1 |
6 | tK | 14 | Ca2+ | 5.4 |
7 | tCl | 15 | CaHCO3+ | 5.4 |
8 | tNO3 | 16 | CaCO3 | 3.0 |
17 | CaSO4 | 3.0 | ||
18 | Mg2+ | 6.2 | ||
19 | MgHCO3+ | 6.2 | ||
20 | MgCO3 | 3.0 | ||
21 | MgSO4 | 3.0 | ||
22 | SO42− | 3.0 |
Total concentrations . | Species concentrations . | |||
---|---|---|---|---|
Index . | Component . | Index . | Component . | Ion size parameter ![]() |
1 | tC | 9 | H+ | 6.8 |
2 | tCa | 10 | OH− | 4.2 |
3 | tMg | 11 | CO2 | 3.0 |
4 | tSO4 | 12 | HCO3− | 1.5 |
5 | tNa | 13 | CO32− | 4.1 |
6 | tK | 14 | Ca2+ | 5.4 |
7 | tCl | 15 | CaHCO3+ | 5.4 |
8 | tNO3 | 16 | CaCO3 | 3.0 |
17 | CaSO4 | 3.0 | ||
18 | Mg2+ | 6.2 | ||
19 | MgHCO3+ | 6.2 | ||
20 | MgCO3 | 3.0 | ||
21 | MgSO4 | 3.0 | ||
22 | SO42− | 3.0 |








Determining equilibrium concentrations (function EqConc)
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)) CO32-, Ca2+, Mg2+ and SO42− are calculated from the given total concentrations. The calculation will be presented taking the example of the carbonate ion CO32−.















Determining pH and SI (functions pH_Alk_Aci and FunctSI)
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−6), 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.
Determining the CCPP (function CCPP_Alk_Aci)
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 CaCO3. Figure 2 shows the flowchart for CCPP_Alk_Aci.
There is no way to directly calculate the CCPP. Therefore, CaCO3 precipitation has to be simulated by iteratively withdrawing virtual amounts of CaCO3 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 CaCO3. For precipitation of X mmol/L of CaCO3, alterations are as follows:
decreases by
, as both CO32− and Ca2+ contribute to the ionic strength with two times their concentration (cf. Equation (23)),
decreases by
, because CO32− contributes to it with two times its concentration (cf. Table 1),
increases by X, because CO32− counteracts CO2-acidity with one time its concentration (cf. Table 1), and
decreases by X, as 1 mol of CaCO3 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−6 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 CaCO3 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.
Main model (model MAIN_CCPP)
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 CaCO3 (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.
RESULTS AND DISCUSSION
Comparison to DIN 38404-10 and a simplified model
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, CaCO3 dissolving. Water 4 is highly dissolving, whereas water 9 is the only one with a positive CaCO3 precipitation potential.
Data of the six example waters and comparison of the pH of calcium carbonate saturation and the CCPP obtained with the developed programme (Progr.) to the values provided in DIN 38404-10 (DIN), and the results of a programme based on a simplified model (Simpl.)
![]() |
![]() |






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 CaCO3 that can be precipitated and underestimate the amount of CaCO3 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 CaCO3. When CaCO3 precipitates, the water's pH decreases, which affects speciation. Consequently, withdrawal of CaCO3 not only results in lower concentrations of Ca2+ and CO32−, but also of other, lower charged species. Therefore, the ionic strength will generally decrease less than four times the precipitated amount of CaCO3. 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.
CONCLUSIONS
This work has shown how to programme calculation of species concentrations, pH, CaCO3 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 CaCO3 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.