Determining freshwater pCO2 based on geochemical calculation and modelling using PHREEQC

Fossil fuel combustion results in rising atmospheric carbon dioxide (CO2), which is known to impact the global climate and the oceans. Latest insights indicate that rising atmospheric CO2 levels also affect CO2 partial pressure (pCO2) in freshwaters, where pCO2 is controlled by a multitude of parameters. However, up to date there is no standardized method, which allows the determination of current and past freshwater pCO2 levels. Ideally methods should incorporate numerous hydrogeochemical and -physical factors to reflect the interplay of all interacting components and their effect on pCO2. We here describe the application of the geochemical program PHREEQC. This freeware serves as an easy method enabling a plausible and comprehensive analysis of pCO2 for field, laboratory, and especially long-term data. We present the use of the different input parameters of a laboratory- and a field long-term monitoring dataset including dissociation constants of carbonic acid measured as total inorganic carbon (TIC) and total CO2 concentration (TCO2) or total alkalinity (TA), together with hydrogeochemical and -physical parameters. Based on current literature and our analyses PHREEQC appears a solid strategy to determine freshwater pCO2 that can moreover be used for long-term datasets.• Comprehensive analysis of pCO2 for field, laboratory, and long-term data.• PHREEQC is not dependent on just one sampling method or parameter scheme.• PHREEQC includes testing the plausibility of a water analysis and enables the assessment of the quality of the laboratory analysis, as well as automatic calculation of all relevant aquatic complexes.


INTRODUCTION/BACKGROUND
With ongoing anthropogenic activities of fossil fuel combustion, CO 2 constantly accumulates in the atmosphere, from where it enters the global carbon cycle [10] . It equilibrates with the ocean water, where it chemically reacts and decreases ocean pH; a phenomenon commonly known as ocean acidification [6 , 11 , 12 , 14 , 17 , 19 , 26] . Similarly, CO 2 accumulation in inland waters was found to be equal or larger than in the ocean [10] . Thus, inland waters have been discussed as hotspots of biogeochemical activity that also act as carbon sinks [10] . Recently, it has been suggested that increased fossil fuel combustion also affects the carbonate equilibria and p CO 2 levels in freshwater systems, which can even lead to decreasing pH ( [10 , 22 , 23 , 25] ; reviewed in [15 , 29] ).
To study this, it is first of all important to chemically describe the freshwater carbonate system. This comprises totally dissolved inorganic carbon content (TIC) that can dissociate into bicarbonate (HCO 3 -), carbonate (CO 3 2-) and aqueously dissolved CO 2 (or H 2 CO 3 * ) [2] . Aqueously dissolved CO 2 consist of two pools, i.e. dissolved free CO 2 and its hydrated form H 2 CO 3 [9] . The concentration of these system components in combination with the calcium and magnesium concentrations, is influenced by the carbonate's equilibria [2 , 27] . All these components regulate pH and determine the buffering capacity [2 , 27] , as well as the entire chemical balance in surface waters [9] . Thus, aqueously dissolved CO 2 is a central variable in aqueous chemistry and is often displayed as CO 2 partial pressure ( p CO 2 ) in a solution [9] . Factors controlling freshwater p CO 2 are atmospheric p CO 2 , the water source (e.g. groundwater, run-off), the residence time of CO 2 in water, the gas transfer velocity (e.g. changes in precipitation rates, wind patterns, etc.) and the underlying waterbody's geology (e.g. limestone or grey weck, etc.) ( [30] ; reviewed in [15] ). This is further complicated by the surrounding soil respiration rate, terrestrial productivity, and land-use but also by biological processes including heterotrophic and autotrophic activities [3 , 4 , 8 , 24] . All these factors directly or indirectly affect the aquatic carbon system and are therefore regarded as p CO 2 -linked-process components (see Table 1 ) rendering p CO 2 analyses in freshwater complicated. For a thorough understanding of dynamics and mechanisms regulating freshwater p CO 2 levels and impacts of climate change linked effects, the analysis of past p CO 2 development is pivotal.
Therefore, a precise determination of current and past freshwater p CO 2 levels, as well as associated p CO 2 -linked-process components calls for a strategy that permits the analysis of field, laboratory and long-term data. Further, an analysis involving not only inorganic but also organic CO 2 -species Table 1 The direct or indirect influence of p CO 2 -linked-process components on freshwater p CO 2 levels. p CO 2 -linkedprocess components are regulating and controlling mechanisms and developments of natural processes. Further contributing factors like climatic and geographical regions, agriculture, as well as hydrogeochemical, and biological processes also influence p CO 2 levels. All these parameters and processes are affected by increased atmospheric p CO 2 by which freshwater p CO 2 is discussed to change. The climate regime has a general influence on the different spheres/processes leading to variation in p CO 2 across freshwater systems and the thereby connected impacts on the p CO 2 . The given p CO 2 -linked-process components can be coded in PHREEQC in the form of keywords or identifiers e.g., the amount of ion species distribution (e.g. Mg 2 + , Ca 2 + , HCO 3 ( [20] ; 2013). With PHREEQC a wide range of equilibrium reactions between water and minerals, ion exchangers, surface complexes, solid solutions, and gases can be simulated [5] . This improves the analysis by reflecting the interplay of all interacting components and their effect on the aquatic p CO 2 .
We here demonstrate the dual use of PHREEQC as a calculation-tool for laboratory datasets based on dissociation constants of carbonic acid. Further, we describe the application of PHREEQC for a field data e.g., stem from long-term monitoring datasets analyzing.

PHREEQC
PHREEQC version 3 is a computer program written in C and C ++ programming and designed to perform a wide variety of aqueous geochemical calculations simulating chemical reactions and transport processes in water [20 , 21] . It is freely available (e.g. https://www.usgs.gov/software/ phreeqc-version-3 ).

METHODOLOGY
With PHREEQC one can evaluate the quality and plausibility of a water sample based on this ion-balance error (i.e. the percent charge-balance error calculated as (100(cationsanions))/(cations + anions) [2] . With this, the positive and negative charges resulting from the single ion charges is determined [21] . As water is electrically neutral, the sum of the positive charges is approximately equal to the sum of the negative charges. Thereby only analyses of water samples with an ion-balance error ≤5% are discussed to be tolerable [2 , 31] . If errors exceed 5% either not all ionic compounds were reflected in the analysis, or individual compound concentrations may be over-or underestimated or individual analysis results are incorrect [2] .

Calculation of pCO 2 and associated complexes
With respect to minerals in a water sample the saturation state ( ) and the activities of free ions in solution can be calculated [2] . Therefore, the solubility product K is compared to the analogue product of the activities derived from water analyses [2] . Further, saturation conditions or the saturation state ( ) may also be calculated as the ratio between IAP (ion activity product of the dissociated chemical species in solution) and K (the solubility product of the mineral) [2] : There is an equilibrium, when = 1, > 1 indicates a supersaturation and < 1 subsaturation. For larger deviations from the equilibrium, a logarithmic scale can be used for the saturation index (SI) [2 , 31] : The SI of different mineral phases, stored in the program's database, is reflected in the PHREEQC output-sheet (see Tab. S7). For dissolved gases the denoted value is given as the common logarithm of partial pressure. Therefore, the calculated p CO 2 is given as the logarithmic value of the saturation index (SI) of CO 2(g) [2 , 31] :

PHREEQC: input of parameters
Program keywords and identifiers can be used for the input of certain hydro-as well as chemical and physical parameters such as pH, temperature, ion species next to the dissociation constants of carbonic acid (e.g. TCO 2 or TA). Parameters obtained from field or long-term data can be entered into the input sheet in following sequence. All given input options are based on Appelo and Postma [2] and Wisotzky et al. [31] . 1 To headline the analysis the keyword ' Title ' is given at the beginning of the input file.
2 Next the keyword ' SOLUTION ' has to be entered directly under the ' Title '. Beneath the keyword ' SOLUTION ', the unit is provided e.g. as ' mg/l ' (Identifier: ' units '). Deviating units, e.g. ' mmol/l ', are written after the corresponding parameter.
3 The parameters temperature, pH, pe and density are provided with the keyword's ' temp ', ' pH ', ' pe ' and ' density ' ( Table 2 ). Pe stands for the Péclet number ( pe = advective transport rate/diffusive transport rate). If density and the pe -value are not available, a default pe ( = 4.0) and a density of 1.0 can be used. The pe of 4 is used in the initial equilibrium calculation for all other redox couples, e.g. of H 2 O/O 2 [2] . 4 Here it is important to enter the measured distribution of carbon in a complete form of carbonic acid species determining p CO 2 and all CO 2 -components [31] . The values for K S4.3 (hydrogen carbonate) and K B8.2 (free carbonic acid) determined by titrations give the total concentration of dissolved CO 2 .This measured distribution of the carbonic acid species can be reflected in the input file in the following ways [31] . The calculations relate to the following information: (a) It is possible to specify CO 2 as carbon without valence: In this case, the program distinguishes the possible oxidation states in C (4) for CO 2 and C (-4) for CH 4 and establishes a redox equilibrium according to the iron specification (reference to the identifier "redox Fe(2)/Fe(3)") [31] . (b) Another possibility is to display the titration results as carbon valence: Here, all carbon in the form of carbonic acid species is given as C (4). By specifying the valence with brackets after the element (C (4)), the speciation of the carbon in methane is explicitly excluded, in contrast to a). Further, if the equilibrium pe-value of the iron redox couple is not in the negative range or it prevails only weakly reducing conditions, the program does not calculate methane concentrations [31] . (c) Alternatively, the hydrogen carbonate determined by titration can be given in the form of alkalinity. In PHREEQC the alkalinity is always stated as CaCO 3 equivalents. The following stoichiometry reaction is decisive for its determination [31] : Only half of the hydrogen carbonate measured in the water comes from the calcite solution, the other half from the carbonic acid. Therefore, the equivalent amount of calcite, based on the acid capacity (K S4.3 ), has to be divided by 2 [31] :

METHOD VALIDATION
We here describe how p CO 2 and all aquatic complexes can be calculated from a water sample of a laboratory dataset and a long-term monitoring dataset with PHREEQC.
The laboratory dataset is based on the dissociation constants of carbonic acid, which is directly measured as total CO 2 concentration (TCO 2 ) and total alkalinity (TA) via endpoint-titration ((K S4.3 ) for TA, acid and base capacity (K S4.3 + K B8.2 ) for TCO 2 ). The dissociation constants TCO 2 and TA are critical input parameters for PHREEQC-analyses as they dependent on water chemistry and buffering capacity of the respective impoundment. In the PHREEQC-input sheet TCO 2 is represented by the keyword C(4)as CO2 (i.e. TCO 2 ) and accordingly as Alkalinity for TA (Table S1; [2 , 31] ).
From a long-term monitoring dataset the annual average of the total inorganic carbon (TIC) can be converted from the total CO 2 concentration. For both data sets, pH and temperature in combination Table 2 A) PHREEQC Input-example in which the parameters can be entered using keywords. Default pe-value ( pe = 4) equating to 230 mV at the given temperature. Total CO 2 given by C (4) as CO2 in mg l −1 = (K S4.3 + K B8.2 (mmol l −1 ) * 44 mg mmol −1 ).
The described alternatives for p CO 2 determination via alkalinity or identify of CO 2 as carbon (without details of the value) are given as Alkalinity (Alkalinity in mg l −1 = (K S4.3 (mmol l −1 ) * 100,089 mg mmol −1 (M CaCO 3 )) / 2) # or C as C (C in mg l −1 = (K B8.2 + K S4.3 (mmol l −1 ) * 12 mg mmol −1 ) [31] . An example sheet for M4 can be found in the supplementary material (S1). B) Example of a model-calculation for the solubility of calcium carbonate and CO 2 in water.

Parameters
Units PHREEQC-Keyword Input-Parameter(to be determined)

A)
Title SOLUTION # or SOLUTION_SPREAD  with a multitude of hydrogeochemical and -physical parameters were included in the p CO 2 -analysis to reflect the interplay of many interacting components and their effect on p CO 2 ( Table 2 ; Tables S1; S5; S6). A problem with long-term data is, that the datasets are often limited in the hydrogeochemical parameters they provide. Hence, we here tested how important the multitude of factors is for the analysis. We therefore performed our analyses with and without the multitude of factors to determine possible deviations between calculated p CO 2 (with an ion-balance error ≤5%) compared to modelled p CO 2 (with an ion-balance error ≤99%). In the following, the abbreviations p CO 2 [TCO 2 ] and p CO 2 [TA] stand for PHREEQC-analyses based on total CO 2 concentration (TCO 2 ) or total alkalinity (TA). TCO 2 measured by endpoint-titration via acid-(K S4.3 ) and base capacity (K B8.2 ) and TA being determined via alkalinity using an endpoint-titration of the acid capacity (K S4.3 ). The symbol ( + ) stands for PHREEQCanalyses including the parameters pH, TCO 2 or TA, temperature, as well as detailed hydrogeochemical and -physical parameters and with a calculated output ion-balance error of ≤5% (listed in Tables 3  and 4 ; Tables S1; S5-S7). The symbol (-) stands for PHREEQC-analyses based on the parameters pH, TCO 2 or TA, temperature, and density only; neglecting hydrogeo-and chemical parameters as well as a calculated output ion-balance error of ≥99%.

Medium
We used a modified version of the artificial M4 medium in all laboratory experiments to ensure a constant water source for the p CO 2 analysis of a water sample (described by the Organization for Economic Cooperation and Development (OECD) in guideline 202 in 2004 [18] ; Table 3 ). This medium contains p CO 2 -linked-process components in the form of the listed ion species (see Table 3 ; S1). Table 3 The concentrations of all elements in the M4 medium of a water sample of the laboratory dataset with modified M4-medium. We freshly prepared 5 L of M4 medium 48 h at the beginning of the experiment, which was split into 4 different experimental conditions of 1 L media. In these we adjusted four different levels of pH (control condition of pH 8.0, p CO 2 -conditions with pH 7.0, 6.7 and 6.4 via CO 2 aeration (Air Liquide-ALIGAL 2)). Elevated p CO 2 -conditions were maintained and monitored with computer-controlled pH probes and solenoid valves (AT-control system, AquaMedic GmbH, Germany) with differences in pH not higher than ± 0.05 units. From each condition we sampled 200 mL for endpoint-titration. Each treatment was independently replicated three times ( n = 12).

pH and temperature
Temperature and pH were determined with a pH electrode and an integrated temperature sensor (LL-Aquatrode plus Pt 10 0 0 F/4 mm; Metrohm AG, Herisau, Switzerland). pH was calibrated using three low ionic strength pH buffers ( We used 0.1 N HCl and 0.1 N NaOH as titrant according to the following equations: where c is the molar concentration of titrant, V 1 the required volume needed to adjust to pH 8.2 or pH 4.3 and V 2 is the volume of the original sample. For higher accuracy we weighed the sample and converted mass to volume using the density of 1 g ml −1 . Results are given in mmol l −1 or mol/m ³ with up to three decimal spaces. TCO 2 was then calculated according to the formula: Similarly, we determined total alkalinity (TA) with 0.1 N HCl for acid capacity (K S4.3 ) via endpointtitration as described above. The concentration of total alkalinity (TA) was then calculated according to the formula: An example of a parameter scheme, analysis results and PHREECQ keywords of a laboratory water sample (M4-medium) is given in the supplement (Table S1). All parameters were calculated according to the composition of modified M4-medium ( Table 3 ) and to the determination-results of acid-(K S4.3 ) and base capacity (K B8.2 ) ( Table 4 ). These parameters were entered via the listed keywords according to the database-sheet in PHREEQC (Tables S1; S7).

Long-term monitoring data
We analyzed p CO 2 over time from data of four freshwater reservoirs in North-Rhine Westphalia. The reservoirs (Henne, Lister, Möhne and Sorpe) are impoundments of approximately 40 m depth and of an age of 61 to 104 years. Long-term monitoring data of hydrogeochemical-andphysical parameters have been monthly recorded since 1970 and data originate from lake surveys (with accredited laboratory for physicochemical water analyses) by the Ruhr Association in Essen (Ruhrverband). All reservoirs are located in the catchment area of the river Ruhr in the western part of Germany and are fed by different tributaries. Further, the area consists of more than 50-60% forest, only about 0.1% moor and about 2-3% are building area. Geographically, the area consists mainly of sandstones, greywacke and claystone (see supplementary information for details).
The conducted p CO 2 analysis based on PHREEQC uses annual averages of the total inorganic carbon (TIC) converted as total CO 2 concentration (TCO 2 ), as well as the annual averages of pH and temperature in combination with hydrogeochemical and -physical parameters from the years 2016-2018 (Table S5).
We determined TCO 2 using the total inorganic carbon (TIC; mg l −1 ) and converted to TCO 2 according to the formula:

Reservoir
Year

POTENTIAL SOURCES OF ERROR
For a given amount of CO 2 in solution, the p CO 2 will change significantly with temperature [9] . Indirect methods using only pH, temperature and TA to calculate p CO 2 are manifold but are highly dependent on the quality of measured parameters [1 , 16] . Small analytical errors may result in large uncertainties in the indirectly calculated TA [1] . In contrast, the quality assessment of a water sample and the laboratory analysis with PHREEQC will help to identify analytic errors of measured parameters.
Especially, in inland waters the choice of TA or TCO 2 used as an input parameter should depend on the respective water chemistry and buffering capacity of the respective impoundment. Therefore, the use of TCO 2 as an input parameter for the dissociation constants of carbonic acid is recommended [1 , 7 , 16] , if data collection or measurement techniques provide TIC/DIC or K B8.2 in combination with K S4.3 (in older datasets also be given as an m-and p-value). These parameters are independent of non-carbonate alkalinity (NCA). Especially in acid, poor buffered inland waters, in which the NCA fraction represents an unknown or even dominant part of alkalinity, it is recommended that estimates of p CO 2 are based on DIC and pH or even direct measurements [1 , 7 , 16] . Here, it is especially important to avoid degassing as otherwise TIC is reduced and CO 2 underestimated [28] . If the type and concentration of the NCA is known, it can be accounted for in the calculation of CO 2 by determining the carbonate contribution as the difference between the total and the non-carbonate alkalinity [28] . Hence, in well-buffered or not acidic and organic rich inland waters and if the NCA fraction is known, also TA is a suitable input parameter. However, input parameters like pH/TA/DIC and water temperature are routinely measured in comprehensive monitoring assays by many environmental agencies [1] . Therefore, these data are readily available to calculate long-term p CO 2 -values helping to detect p CO 2 development processes in the past.

FINAL REMARKS
Here PHREEQC ( + ) is recommended as a calculation method for the standardized analysis of p CO 2 levels giving realistic p CO 2 -calculations of current and past field-, laboratory-and especially long-term data. However, in many long-term datasets there are only a few or even no additional parameters that only allows a modelling of p CO 2 and will result in analyses with an ion-balance error of ≥ 99%. Nevertheless, we here find, that modelled p CO 2 can still be used as an estimation close to the actual p CO 2 level.
Each freshwater body is unique, due to chemical compositions of the underlying geology, climate regimes and/or freshwater chemistry and therefore the p CO 2 and pH relationship is less straight forward compared to the marine system, rendering analyses more difficult. For a global comprehension of freshwater p CO 2 , future field and laboratory experiments are needed to further describe p CO 2 dependent freshwater acidification also for different geographic locations. Foremost, this requires a standardized methodology with which current and past p CO 2 can be determined. In future research PHREEQC could be used to analyze the development of current p CO 2 in freshwater from different sources like rivers, lakes and reservoirs capitalizing on long term dataset from global databases (e.g. GloRiCh, Global Reservoir and Dam Database (GRanD), Waterbase (EEA's databases), waterdatabase) or environmental agencies. Further, the advantage of PHREEQC is that with just one method p CO 2 in very different sample types can be analyzed. Here, PHREEQC can help to provide a deeper understanding of p CO 2 regulating and influencing processes predicting future effects of climate change.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.