Exploring uncertainty in thermodynamic modeling of the lead carbonate aqueous system

Lead is a pervasive toxic metal that can exist in the home, air, soil, or water. Industries focused on the removal or control of lead include wastewater treatment, battery recycling, and soil remediation. Furthermore, lead-bearing plumbing materials in drinking water distribution systems are common. The release of lead from corrosion scales that form inside these plumbing materials can cause elevated lead concentrations in drinking water at the tap. While many municipalities have programs for the replacement of lead-bearing plumbing materials, these programs are costly and progress slowly (Sandvig et al. 2008). As such corrosion control strategies are often implemented to reduce lead concentrations in drinking water. Water chemistry is generally the most accessible parameter for water utility operators to adjust to control corrosion and thus reduce lead concentrations. For instance, parameters such as alkalinity (or dissolved inorganic carbon, DIC) (Schock 1980; Marani et al. 1995; Tam & Elefsiniotis 2009; Wang et al. 2012); pH (Tam & Elefsiniotis 2009; Kim et al. 2011; Wang et al. 2012); additives such as orthophosphate (Gouider et al. 2009; Tam & Elefsiniotis 2009; Ng et al. 2012); and disinfectant type (i.e., free chlorine and chloramines) (Edwards & Dudi 2004; Vasquez 2005; Switzer et al. 2006; Lin & Valentine 2008; Lin & Valentine 2009) can be adjusted to stabilize lead-bearing corrosion scales. As conditions vary between distribution systems, there is a need to be able to predict a priori the effect water chemistry adjustments may have on lead corrosion scale stability and thus lead concentrations.

The composition of the corrosion scale is important for determining the effectiveness of a corrosion control strategy (Docherty & Kariuki 2019). The main groups of phases in lead-bearing corrosion scales are lead (II) carbonate and lead (IV) oxide solids (Wang et al. 2012). The major lead (II) carbonates are hydrocerussite (Pb₃(CO₃)₂(OH)₂) and cerussite (PbCO₃), while the major lead (IV) oxides are plattnerite (PbO₂) and scrutinyite (PbO₂) (Wang et al. 2012). The characterizations of experimental pipe loops and exhumed lead pipes have shown that lead (II) carbonates are often the dominant solid phase in distribution systems with higher DIC (Edwards et al. 1999; Liu et al. 2008; Noel et al. 2014).

Numerous studies over the last 50 years have used thermodynamic models to describe lead carbonate aqueous–solid phase interactions for varying water chemistry conditions. Details of models previously used are provided in Supplementary Material, Table S1. These thermodynamic models have provided valuable insight into the relationship between water chemistry, solid phases, and lead solubility (Hem & Durum 1973; Patterson et al. 1977; Schock 1980; Marani et al. 1995; Noel & Giammar 2008; Tam & Elefsiniotis 2009; Noel et al. 2014; Tully et al. 2019). A common theme among these studies is uncertainty in the lead solid phases included in the model (i.e., phases allowed to precipitate) and the equilibrium constants (logK values) used for the lead solid phases and aqueous complexes. There is a large variation in the values for equilibrium constants reported in the literature, and this can result in lead solubility and solid phase predictions that differ over several orders of magnitude between thermodynamic models, and between model predictions and experimental or field data (Noel & Giammar 2008; Noel et al. 2014). Equilibrium constants are typically determined by measuring species of interest under different experimental conditions (Vanbriesen et al. 2009). Equilibrium constant values reported may vary due to multiple reasons including the specific conditions of the experiment (temperature, pressure, ionic strength, and solution species present), the analytical methods used to measure the species concentrations (with methods improving over time), or variations in the morphology of the solid phases.

Based on prior thermodynamic models and laboratory experiments (see Supplementary Material, Tables S1 and S2), there are two dominant solid phases (hydrocerussite and cerussite) and up to nine lead aqueous species (⁠⁠, Pb(OH)₂, Pb(OH)₃⁻, Pb₂OH⁺³, ⁠, ⁠, PbCO₃^o, PbHCO₃⁺, and PbOH⁺) in the basic lead carbonate system. Based on thermodynamic constraints, hydrocerussite is the most stable lead (II) carbonate when the pH is above 7–7.5, whereas cerussite is stable under more acidic conditions (Schock 1980; Schock et al. 1980). A range in equilibrium constants has been reported for hydrocerussite and cerussite. For instance, the hydrocerussite equilibrium constant (logK) has been reported to range from −16.25 to −19.04 (Supplementary Material, Tables S1 and S2). While Noel et al. (2014) and Tully et al. (2019) recently illustrated the sensitivity of lead solubility predictions to the hydrocerussite equilibrium constant, predictions also depend on the equilibrium constants for lead aqueous species. Equilibrium constants for these species also have a range of reported values. There is a need to understand the way in which uncertainty in equilibrium constants, including those for aqueous species, propagates through the thermodynamic model to guide the selection of equilibrium constant values, understand uncertainty in predictions, and thus to develop more robust thermodynamic lead solubility models.

Uncertainty in thermodynamic models caused by variation in equilibrium constant values has been explored previously for other chemical systems. For instance, Vanbriesen et al. (2009) used a Markov chain Monte Carlo to regress 27 thermodynamic constants for an EDTA system simultaneously. More recently, Guo et al. (2019) used a Monte Carlo with Latin Hypercube sampling to propagate equilibrium constant uncertainty through thermodynamic simulations of the mercury system. Both studies found that slight variations in equilibrium constants can result in significant changes to the output species concentrations spanning multiple orders of magnitude. Furthermore, Guo et al. (2019) highlighted that there is a significant lack of thermodynamic parameter uncertainty quantification for environmentally important metals such as As, Pb, Cd, and Hg.

While Monte Carlo methods are popular, they require probability distributions for input parameters (Novoselov et al. 2015). If probability distributions for input parameter values are not available, as is the case for the equilibrium constants for the lead carbonate system, options to evaluate uncertainty are more limited. Modifying one factor (i.e., one equilibrium constant value) at a time can be cumbersome and miss critical interaction effects between factors. The factorial analysis is a sensitivity analysis method useful for assessing many factors without requiring large computational resources (Liu et al. 2020). A full factorial analysis involves simulating the permutation of every combination of factors at every level (high or low value) investigated. One of the most important aspects of factorial analysis is the ability to analyze interactions between the factors (Mee 2009). For example, Liu et al. (2020) showed important interaction effects between factors (e.g., between snowfall temperature and threshold temperature for snowmelt) in a model used to simulate average streamflow. In a thermodynamic model, these factors would be the equilibrium constants for the various species included in a model for simulating total dissolved lead.

The main objective of this study is to evaluate how uncertainties in equilibrium constants for lead solid phases and aqueous complexes propagate through a thermodynamic model used to simulate the lead carbonate aqueous system. In particular, the study aims to (i) identify species' equilibrium constants that contribute most to model uncertainty; (ii) determine how uncertainty contributed by each species equilibrium constant varies with changing pH and DIC conditions; and (iii) provide recommendations for future thermodynamic modeling of this system. The lead carbonate system considering two solid phases (cerussite and hydrocerussite) and nine aqueous lead complexes is simulated using a full factorial statistical approach that considers the range of equilibrium constant values previously reported. Batch laboratory lead solubility experiments are also conducted with results compared with output from the factorial analysis. While the lead carbonate system simulated in this study is simplified relative to the water chemistry and solid corrosion scales found in drinking water distribution systems, analyzing this simpler system provides important fundamental insight into the uncertainty inherent in thermodynamic simulations of the aqueous lead system. The study findings can be applied to guide the future application of thermodynamic models including consideration of uncertainty inherent in these models.

A full factorial statistical analysis was performed to analyze the way in which uncertainty in equilibrium constants (logK values) for all lead species in the lead carbonate system propagates through the thermodynamic model. A custom Python (v.3) wrapper code was used to iterate through factorial combinations of logK values that were used as parameter inputs to a thermodynamic model implemented in PHREEQC (v1.4.2). This code has been uploaded to a public repository (see Data Availability Statement for URL). The PHREEQC model was a batch simulation with the lead carbonate solid-phase hydrocerussite initially present in excess in a buffered solution with ionic strength 0.01 M that was closed to the atmosphere at room temperature (20 °C). A second lead carbonate, cerussite, was not initially present but was allowed to precipitate if the saturation index is >0. Nine aqueous lead species were considered: ⁠, Pb(OH)₂, Pb(OH)₃⁻, Pb₂OH⁺³, ⁠, ⁠, PbCO₃^o, PbHCO₃⁺, and PbOH⁺. These were chosen based on prior modeling and experimental studies (Schock 1980; Marani et al. 1995; Xie et al. 2010a, 2010b; Noel et al. 2014). Dissolution and formation equations for the lead aqueous species and solid phases are provided in Table 1. The MINTEQ.V4 database (Allison et al. 1991) was used to specify the equilibrium constants for non-lead aqueous and solid species (e.g., H₂CO₃, HCO₃⁻, and ⁠) that were considered in the model. Simulations were conducted with pH ranging from 4 to 11 and DIC ranging from 0 to 250 mg C/L. These pH and DIC conditions were coded into the PHREEQC model via the python wrapper.

The custom python wrapper substituted logK values into the PHREEQC batch simulation for each solid and aqueous lead species. This substitution was based on a target logK table (Table 1 and Supplementary Material, Table S3) which identified the levels (e.g., minimum and maximum) for each factor to be included in the factorial analysis. The levels used for each logK value were based on an extensive literature review including Lothenbach et al. (1999) who previously compiled relevant logK values for the lead system that were reported in studies published between 1922 and 1998. Primary references for logK values obtained from Lothenbach et al. (1999) are provided in Supplementary Material, Tables S1 and S2. Once all combinations of logK values were then run through the PHREEQC simulation, the python code collated the model output for analysis. An initial factorial analysis was run with 11 factors 2 levels (2¹¹), meaning that each of the 11 factors (aqueous or solid phase logK values) investigated were varied across 2 levels: a high level and a low level. For this factorial analysis, a level was either the minimum or maximum logK value previously reported in the literature for that factor (Table 1).

A second factorial analysis was run with 5 factors 12 levels (12⁵). This factorial analysis only included the five factors that were identified in the 2¹¹ analysis to be associated with the greatest uncertainty in the model output. Additional levels that fell between the minimum and maximum values were included in this factorial analysis to further understand the sensitivity of each factor (Supplementary Material, Table S3). In this factorial analysis, the logK values in the MINTEQ. V4 database (Allison et al. 1991) were also set as a level for each factor. Finally, the results of the 12⁵ factorial analysis were compared with batch lead solubility experimental results to evaluate whether the experimental data could be simulated when considering all logK combinations for these five factors that contribute to the greatest uncertainty in the model output. In doing this, recommendations are provided for selecting logK values for these five factors and for considering uncertainty in future thermodynamic simulations of the lead carbonate system.

The model output for the 2¹¹ factorial was also analyzed to examine interaction effects. It is important to understand interaction effects because they show whether the selection of one factor (logK value for one species) may change the uncertainty contribution of another factor (logK value for a different species). If interaction effects are high between different factors, this indicates that it is not possible to select the logK value for an individual species without considering the way in which the value depends on the logK value selected for another species. This would mean that the selection of appropriate logK values for all species in an aqueous system may need to be done simultaneously.

Batch hydrocerussite dissolution experiments were conducted across drinking water conditions (pH 8–10; DIC 20, 100, and 200 mg C/L) to determine whether the thermodynamic model was able to simulate total lead solubility from hydrocerussite dissolution when considering uncertainty in the logK values. The matching of the experimental data with the 12⁵ factorial model output was used to determine appropriate logK values (or combinations of values) for the aqueous and solid-phase species to be used in the model. It should be noted that while the logK values selected based on the matching are appropriate for the solid phases used in our experiments, it is possible that logK values for solid phases may vary with crystallite size and other aging phenomenon.

Hydrocerussite (Pb₃(CO₃)₂(OH)₂; Sigma-Aldrich) was used for all experiments. XRD analysis and Rietveld refinement indicated that the solid phase was dominant hydrocerussite. Sodium bicarbonate (NaHCO₃; Sigma-Aldrich, >99.5%) and sodium nitrate (NaNO₃; Sigma-Aldrich, >99.0%) were combined with mega-pure water to prepare a stock solution at the target DIC and at 0.01 M ionic strength, respectively. The target pH was maintained with organic buffers: EPPS (Alfa Aesar, 0.2 M, pH 8), CHES (Alfa Aesar, 0.5 M, pH 9), CAPSO (Alfa Aesar, 0.2 M, pH 9.5), and MOPSO (Alfa Aesar, 0.2 M, pH 7). Sodium hydroxide (NaOH; Sigma-Aldrich) and concentrated nitric acid (HNO₃; EMD Milipore, Omnitrace, 67–70% v/v condensed) were used to adjust the stock solution to the target pH (±0.05 pH units) immediately before each experiment.

All experiments were performed as closed experiments at room temperature (20 °C). The influence of the initial amount of hydrocerussite present (solid loading) on the time to equilibrium was first evaluated by conducting experiments with varying solid loadings (10–1,000 mg/L; Supplementary Material, Figure S1) at pH 8 and DIC 20 mg C/L. Equilibrium was determined to have been reached when the measured dissolved lead concentrations were within 10% of each other for at least three samples collected at consecutive times. Based on these experiments, a solid loading of 175 mg/L was used for all subsequent experiments with a minimum of 4 days required until equilibrium was reached.

Experiments were run at conditions typical for drinking water systems and for which hydrocerussite is the dominant solid phase (pH 8–10) and at varying DIC levels (20, 100, and 200 mg C/L; Supplementary Material, Table S4). Complications due to possible cerussite formation can occur at pH lower than 7.5 (Schock 1980). Therefore, while some drinking water systems operate at lower pH, all experiments were run at pH 8 or higher to avoid complications of cerussite/hydrocerussite dominance within the experiments. Organic buffers were used to maintain target pH values, since a positive pH-driven feedback loop can occur in unbuffered solutions as hydrocerussite dissolution increases as pH increases. The buffers used have been used in previous hydrocerussite dissolution experiments (Noel et al. 2014). Additional experiments (not presented) were conducted which indicated that organic buffers do not alter the dissolved equilibrium lead concentrations at the range of pH values studied.

For all experiments, powdered hydrocerussite was first weighed into 150 mL polypropylene bottles. Each bottle was then half-filled with a pre-prepared stock solution and agitated to evenly disperse the hydrocerussite powder. After agitation, the bottles were filled completely and placed on a shaker table. A minimum of five sample bottles were run for each experimental data point (Supplementary Material, Table S4). While the experiments were run under closed conditions, the stock solution was not prepared under closed conditions and so alkalinity was measured to verify target DIC values. Interference from atmospheric CO₂ was not observed with all DIC measurements at expected target concentrations.

All samples were taken from individual bottles that were sacrificed after sampling. Samples were immediately filtered through a 0.2 μm polyethersulfone filter. It is possible that the selection of this filter size (0.2 versus 0.45 μm) may have influenced the dissolved equilibrium concentrations and thus the factorial analyses matching results, but the influence of filter selection was not examined in this study. Samples for total dissolved lead analysis were acidified immediately using 2% HNO₃ and analyzed using graphite furnace atomic absorption spectroscopy (GFAAS, Agilent Technologies; limit of detection for lead = 0.001 mg/L). Standard lead drying, ashing, and volatilization procedures for lead were followed without a matrix modifier. Experimental blanks and laboratory blanks and spikes were included in analysis for quality control and assurance.

The total dissolved lead concentrations for the different experimental conditions were compared with the 12⁵ factorial model output to determine set(s) of logK values that could be used to simulate the experimental data within a given tolerance. The tolerance was set to five times the standard deviation for each experimental data point (Supplementary Material, Table S4). When a factorial model output of total dissolved lead was within the tolerance of an experimental point, the grouping of logK values used for that factorial PHREEQC simulation were considered a ‘match’. Common sets of logK values that ‘matched’ all nine experimental data points were identified and assembled into groups for evaluation (Table 2).

The main effect of a factor (species logK value) indicates the impact of this factor on the dependent variable (total dissolved lead). A large main effect indicates that changing the factor from a low value to a high value has a large effect on the simulated total dissolved lead concentration. The model uncertainty associated with a factor with a large main effect is high. The calculated average main effects for all factors considering the entire condition space simulated (pH 4–11, DIC 0–250 mg C/L) are shown in Figure 1(a) with the data provided in Supplementary Material, Table S5A. Hydrocerussite (22.3 mg Pb/L), (12.4 mg Pb/L), PbOH⁺ (7.6 mg Pb/L), cerussite (4.1 mg Pb/L), and PbCO₃^o (1.1 mg Pb/L) were found to have highest average (across the entire condition space) main effects. The total sum of average main effects considering all factors is 48.7 mg Pb/L – this value can be helpful in determining the relative importance of each factor on the model output (total dissolved lead). In comparing this main effect to the magnitude of experimental lead concentrations (<1 mg/L), the importance of logK value selection is evident, as the main effect of a selected logK can have a large impact on the simulation results when compared to experimental concentrations.

Hydrocerussite accounts for 46% of the average main effects, indicating that almost half of the uncertainty of the model arises from uncertainty in its logK value. represents 25% of the average main effects, with PbOH⁺ (16%), cerussite (8%), and PbCO₃^o (2%) accounting for the remaining approximately 25% of the average main effects. All other factors had average main effects of less than 1 mg Pb/L and account for only approximately 3% of the total average main effects, indicating that uncertainty associated with these factors is minimal. Based on average main effect analysis, there are two primary factors: hydrocerussite and and three secondary factors: PbOH⁺, cerussite, and PbCO₃^o that contribute to uncertainty in the total dissolved lead predictions. These five factors (hydrocerussite, ⁠, PbOH⁺, cerussite, and PbCO₃^o) were further investigated by running a 12⁵ factorial analysis.

While average main effects indicate the influence of changing one factor across the entire condition space on the dependent variable (total dissolved lead), it is also helpful to examine the maximum main effects for each factor and at what conditions they occur (Figure 1(b), Supplementary Material, Table S5B). The maximum main effect is the highest main effect for each factor considering each condition (pH 4–11, DIC 0–250 mg C/L). The largest maximum main effects were for (146 mg Pb/L, pH 9.5, DIC 250 mg C/L), hydrocerussite (95.5 mg Pb/L, pH 9.5, DIC 250 mg C/L), cerussite (86.7 mg Pb/L, pH 4, DIC 250 mg C/L), and PbOH⁺ (28.8 mg Pb/L, pH 4, DIC 0 mg C/L). This indicates that model uncertainty due to the logK value was the largest and can change the simulated total dissolved lead by up to 146 mg Pb/L at a pH value of 9.5 and a DIC value of 250 mg C/L. The maximum main effects all occur near the limits considered in the study, typically at either pH 4 or pH 10.5. It is therefore recommended that greater care be taken in selecting equilibrium values for all factors when simulating these extreme conditions. All other factors studied had maximum main effects less than 6 mg Pb/L with these maximum main effects again occurring at extreme pH values (4 or 10.5).

The main effects across the condition space for the five factors identified to have the largest average main effects are shown in Figure 2 (see Supplementary Material, Figure S2 for other factors). Note that the heat maps have an uneven contouring scale where green is less than 10 mg Pb/L. Typical drinking water conditions (between 20 and 200 mg C/L and pH 7 and 10) are also indicated. From Figure 2(a), it can be seen that uncertainty in the hydrocerussite logK results in three main areas of high model uncertainty: low pH (<5.5); mid-high pH (8–10) and high DIC (100–250 mg C/L); and extremely high pH (11). At these conditions changing the hydrocerussite logK value has a large impact on the total dissolved lead predictions. Furthermore, the results suggest that any model simulating drinking water conditions will have high uncertainty associated with the selection of the hydrocerussite logK value. Cerussite has a similar main effect heat map to hydrocerussite with three regions of high uncertainty: low pH (<5) with mid to high DIC (>50 mg C/L); mid-high pH (9–10) with mid to high DIC (>100 mg C/L); and extremely high pH (11; Figure 2(b))). For typical drinking water conditions (pH>7 and DIC >0 mg C/L), and there is moderate uncertainty associated with the selection of the cerussite logK value. has one area that has high main effects at mid to high pH and mid to high DIC – this overlaps with typical drinking water conditions so care needs to be taken in selecting this logK value (Figure 2(c)). Finally, PbOH⁺ and PbCO₃^o have low main effects when considering typical drinking water conditions with both species having higher main effects at lower pH (pH <8) or at extremely high pH (11; Figure 2(d) and 2(e)).

Interaction effects indicate whether a factor is dependent upon another factor. This is important to consider because if there is a large interaction between factors, both factors need to be simultaneously considered/adjusted when selecting their values as their impact on the model output is interconnected. Prior thermodynamic modeling of the lead carbonate system has only examined the effects of modifying a single logK at a time to determine its impact on the model output, thereby neglecting interaction effects (e.g., Noel et al. 2014; Tully et al. 2019).

Large interaction effects observed from the factorial analyses are shown in Figure 3. In interpreting these figures, parallel lines indicate no interaction between factors, while divergent slopes indicate high interaction (results for all interaction effects are provided in Supplementary Material, Table S6). As can be seen in Figure 3(a)–3(c), cerussite has three major factors that it interacts with: ⁠, PbOH⁺, and PbCO₃^o. As such the cerussite logK value needs to be considered in assessing the effect of the logK values for these other species. For instance, as the logK for cerussite increases, the model uncertainty associated with the logK values of (Figure 3(a)) and PbCO₃^o (Figure 3(c)) increases. Alternatively, as the logK for cerussite increases, the model uncertainty associated with the PbOH⁺ logK value decreases (Figure 3(b)). These interaction effects mean that it is important to determine the logK values for cerussite, ⁠, PbOH⁺, and PbCO₃^o simultaneously. Cerussite only has minor interactions with other species including Pb(OH)₂, ⁠, ⁠, and PbHCO₃⁻ (Supplementary Material, Table S6). Finally, hydrocerussite has two major interactions: (Figure 3(d)) and PbOH⁺ (Figure 3(e)) with the uncertainty associated with both these species logK values increasing as the hydrocerussite logK increases. A minor interaction between hydrocerussite and PbCO₃^o is shown in Figure 3(f) for comparison. This indicates that the logK values used for hydrocerussite, ⁠, and PbOH⁺ should be determined simultaneously.

The results for the lead solubility experiments at various pH values (8–10) and DIC values (20, 100, and 200 mg C/L) are shown in Figure 4. All experimental points at pH 8 have a similar total dissolved lead concentration (∼0.05 mg/L) regardless of the DIC with the lead concentrations diverging as pH increases. For the conditions examined, the total dissolved lead is the highest for pH 10 and DIC 200 mg C/L (0.62 mg/L).

Before comparing the experimental data to the 12⁵ factorial results, the contribution of each of the factors considered in the 2¹¹ factorial analysis to the total main effects was examined to confirm that the five factors included in the 12⁵ factorial analysis accounted for the majority of the model uncertainty at the experimental conditions. Of the factors considered in the 2¹¹ factorial analysis but not included in the 12⁵ factorial analysis, only Pb(OH)₂ was found to contribute to uncertainty (<4%) at low DICs (see Supplementary Material, Figure S3). Twenty-one combinations of logK values out of the 248,832 combinations considered in the 12⁵ factorial analysis were found to fit all experimental data points with a tolerance of five times the standard deviation between replicate experimental samples. These 21 combinations of logK values that match experimental data can be grouped according to the logK values of ⁠, hydrocerussite, and cerussite to form six general groups (Table 2). The uncertainties in the logK values for PbOH⁺ and PbCO₃^o have less effect on the model output relative to the logK values for hydrocerussite, cerussite, and ⁠. As the logK value for PbOH⁺ did not considerably affect the dissolved lead concentrations, a range of values could be used.

Two parameters, the average difference and the maximum difference between the experimental values and the simulated values, were calculated to quantitatively compare the goodness of the fit of the model for each of the general sets of logK values (Table 2). The average difference considering all the experimental conditions ranged from 0.021 mg/L (Set 3) to −0.0686 mg/L (Set 6), whereas the maximum difference ranged from 0.101 mg/L (Set 2) to −0.432 mg/L (Set 4). The logK values in Sets 1a and 1b provided the closest match with the experimental data considering the calculated average and maximum differences (Table 2). While the matching suggests that any of the combinations of logK values shown in Table 2 could be used in future simulations, in general, it is recommended that the logK values associated with Set 1 (Table 2; Figure 4) are used as it was the most common result from matching. In contrast to prior studies that have been unable to simulate experimental results within the tolerance used in this current study (e.g. Noel et al. 2014), here this was possible by understanding the uncertainty and the interactions between the logK values for the multiple solid and aqueous phase species in the lead carbonate system. This highlights the need to consider the uncertainty of all logK values and the interactions between these values when using thermodynamic models to simulate the lead carbonate aqueous system.

This study has shown how uncertainty in equilibrium constant values for solid and aqueous phase species propagates through thermodynamic models of the lead carbonate aqueous system. It was determined that the logK values for the following species contributed most to uncertainty: hydrocerussite, ⁠, cerussite, PbOH⁺, and PbCO₃^o. The model output was shown to be highly sensitive to the hydrocerussite logK value with average main effects, indicating that this value accounts for 46% of the uncertainty in the thermodynamic lead solubility predictions for pH 4–11 and DIC 0–250 mg C/L. When simulating typical drinking water conditions (between pH 7 and 10 and DIC 20 and 200 mg C/L), care must also be taken in selecting the equilibrium constant for ⁠, as it showed the most effect under these conditions. The logK value was associated with the second-highest average uncertainty (up to 25%) across the pH and DIC condition space considered, and also had the highest maximum main effect. Overall, maximum main effects indicated that model uncertainty due to logK selection is typically greatest at extreme pH values (pH <5 and >10). Significant interactions effects were also observed between the logK values for cerussite, ⁠, PbOH⁺, and PbCO₃^o and between the logK values for hydrocerussite and and PbOH⁺. This indicates that logK values used for these species should be selected simultaneously for thermodynamic modeling of the lead carbonate system. Combinations of logK values that provided a good match with experimental data were determined with the average difference between the simulated and experimental lead concentrations calculated to be 0.031 mg/L for the best logK value combination (Set 1a/1b).

Overall, this study provides valuable insight into the uncertainty present in thermodynamic models of the lead carbonate system. While the simulated lead carbonate system is simplified relative to the water chemistry in drinking water distribution systems, analyzing this simpler system provides important insight into the uncertainty inherent in thermodynamic simulations. Lead chemistry is more complex in drinking water distribution systems due to reaction kinetics, temperature changes, and the presence of species including chlorine, natural organic matter, and calcium. It is recommended that this uncertainty analysis be extended to other species including lead phosphate (inorganic) species and that the effects of temperature are explored. Additionally, it is recommended that lead solubility experiments that use well-characterized lead corrosion scales from real distributions systems be combined with similar thermodynamic model uncertainty analysis to further evaluate the sets of recommended logK values derived in this study.

All relevant data are included in the paper or its Supplementary Information. A compilation of logK values from literature sources, tabulated experimental results, additional simulation results, and a summary of all interaction effects are provided in the Supplementary Materials. The python code is used to implement the factorial analysis and matching process, and all main and interaction effects plots are available in a data repository: https://github.com/UWO-RESTORE/Kushnir-2021.