Factors Governing Site and Charge Density of Dissolved Natural Organic Matter

1. Introduction

The concentrations of dissolved natural organic matter (DNOM) have increased in most northern freshwaters. This is attributed to several factors including the reduced impact of acid rain [1,2], a warmer and wetter climate [3,4], and growing biomass [5,6]. Many acid-sensitive surface waters recovering from acid deposition show a return from anthropogenic acidification to natural acidity, where increased anionic charge from DNOM contributes considerably to the total anionic charge. This decline in acid rain, alongside the other contributing factors, has heightened the relative and absolute significance of the organic charge stemming from the corresponding bases of weak acid functional groups present on DNOM.

Organic charge density (CD), which in this study refers to the concentration of calculated organic charge (Org.^-) normalized to the concentration of total organic carbon (TOC), is a conceptually important factor in understanding the solubility of DNOM and its function as a counter ion in the leaching of base cations, and thus recovery from acid rain, as well as its complexation and detoxification of heavy meals and aluminum (Al). Additionally, organic charge plays an important role in the quality control of major anions and cations data, as it contributes significantly to the charge balance control, especially in dystrophic waters.

A modelled estimate of the organic charge concentration (OAN^-) can be derived through conceptually formulated and empirically fitted models based on TOC, pH, and an estimate of the density of weak acid functional groups on the DNOM (SD). However, a challenge arises as this SD is hypothesized to differ in space and vary with time, influenced by changes in the quality of DNOM as well as trends in the concentration of organically complexed aluminum. Consequently, there is a pressing need to elucidate the factors governing SD and CD and to gain a deeper understanding of their influence on these variables in space and time.

2. Theory

2.1. Charge and Site Densities

Charge Density (CD) is determined by normalizing Org.^- to the TOC concentration (i.e., $\frac{{Org.}^{·}}{TOC}$). The Org.^- is determined through ion balance calculations involving the molar equivalent charge of major cations and anions (i.e., calcium (Ca²⁺), magnesium (Mg²⁺), sodium (Na⁺), potassium (K⁺), ammonium (NH₄⁺), protons (H⁺), labile aluminum (LAlⁿ⁺), sulphate (SO₄^2-), nitrate (NO₃^-), chloride (Cl^-), and bicarbonate (HCO₃^-)). In cases where concentrations of K⁺, NH₄⁺, LAlⁿ⁺ and HCO₃^- are absent and negligible, they may be excluded. However, if the sample pH exceeds 5.5, a bicarbonate measure based on alkalinity readings must be factored into the calculation.

Brakke et al. [7], analyzing ion balance in water chemistry data from 94 dystrophic lakes across two regions of Norway (i.e., Gulsvik and Nord-Trøndelag), proposed the empirically fitted linear Eqn. 1 for estimating organic charge.

Org.^- = CD ∙ TOC – 0.96

(1)

They stated the CD value as 4.06 µEq/mg C. This straightforward equation accounted for 55% of the variance in Org.^-. Similarly, Henriksen and Seip [8] reported a CD of 5.5 µEq/mg C based on extensive observations in lakes across Norway and Scotland in the late seventies, while Kortelainen [9] found a median CD of 7.5 µEq/mg C studying 955 Finnish lakes sampled in 1987.

The Site Density (SD) is employed in conceptual organic acid models to characterize the relative abundance of weak organic acid functional groups (i.e., mainly carboxylic and phenolic acids) within the DNOM. SD can be determined through advanced measurements of the content of these acids in DNOM. Alternatively, as undertaken in this study, it can be empirically calibrated within the models by optimizing the model simulation of OAN^- relative to calculated Org.^- across extensive datasets.

Based on Gran titration of ion-exchanged samples Hruška et al. [13] observed SD to be 10.2 µEq/mg C for Swedish sites and 8.8 µEq/mg C for the more acidified Czech sites. Oliver et al. [12] proposed an SD of 10.5 µEq/mg C (± 1.0). This was derived from a weighted carboxylic content in humic and fulvic DNOM fractions in a limited number of XAD-8 isolated Canadian samples. Perdue and Richie [14], in a literature review of direct titrations of carboxylic and phenolic acid content in DNOM fractions, found lower mean densities of 4.0 and 1.8 µEq/mg C for carboxylic and phenolic acid groups in humic acids, respectively, compared to 5.1 and 1.7 µEq/mg C for fulvic acids, respectively. Given DNOM’s blend of humic and fulvic acids, the concentration of carboxylic and phenolic acid of the DNOM was measured at 5.2 and 1.9 µEq/mg C, respectively.

Within the fulvic DNOM fraction Perdue and Richie [14] suggested a higher SD of weak acid functional groups, such as carboxylic, and phenolic acids, compared to the humic fraction. Additionally, they documented a greater level of aromaticity in the humic fraction compared to the fulvic fraction. Aromatic compounds, rich in conjugated double bonds, absorb UV light and wavelengths (λ) into the visible (VIS) light spectrum. λ of 254 and 410nm are thus widely used as proxies for DNOM, and their specific absorbencies (sUVa and sVISa) are positively correlated to DNOM aromaticity and, consequently, to the ratio of humic to fulvic fractions in the DNOM ([15] and references therein). Based on this conceptual framework, one could hypothesize that the SD of DNOM decreases with increasing sUVa and sVISa.

In summary, SD and CD of DNOM exhibit significant spatial and temporal variability and their estimates depend on the methodology employed for their determination.

2.2. Organic Charge Models

Assuming a triprotic model, Lydersen et al. [10] proposed the Eqn. 2 for OAN^- by simply assuming one-third of the SD as strong acid anions.

$${OAN}^{·}=\frac{SD·TOC}{3}$$

(2)

This equation has been utilized in critical load assessments [11] to improve the ecological relevance of the acid-neutralizing capacity (ANC) using organic acid-adjusted ANC (ANC_OAA = ANC – OAN^-). Still, the actual mole fraction (α) of charged functional groups hinges on the protonation level of the organic acids and their complexation with aluminum and partly iron (Fe). Consequently, the CD is governed by the SD and influenced by factors such as pH and metal complexation. SD is thus higher than or equal to CD due to the presence of protonated sites or complexes with Al and Fe.

Prominent models accounting for the effect of pH variations on CD include the Oliver model [12] and the Hruška model [13]. These conceptual frameworks simulate the concentration of organic charge (OAN^-) based on TOC, pH, and a selected SD. The Olivers model [10] is given in Eqn. 3.

$${OAN}^{·}=TOC·SD·\frac{{10}^{·\left(0.96+0.9·pH·0.039·{pH}^2\right)}}{{10}^{·\left(0.96+0.9·pH·0.039·{pH}^2\right)}+{10}^{·pH}}$$

(3)

Here the fraction expression denotes the mole fraction (α) for the overall protolysation ratio (i.e., CD/SD). Spatial disparities and temporal changes in SD may be accommodated by adjusting this parameter within the models.

The Hruška model is based on the definitions of mole fractions (α) of the protonated acids considering the organic acids on the DNOM as triprotic (Eqns. 4-7).

$$\alpha 0=\frac{{\left({10}^{·pH}\right)}^3}{{\left({10}^{·pH}\right)}^3+{10}^{·pKa1}·{10}^{·pKa2}·{10}^{·pH}+{10}^{·pKa1}·{10}^{·pKa2}·{10}^{·pKa3}}$$

(4)

$$\alpha 1=\frac{{10}^{·pKa1}·\alpha 0}{{10}^{·pH}}$$

(5)

$$\alpha 2=\frac{{10}^{·pKa1}·{10}^{·pKa2}·\alpha 0}{{\left({10}^{·pH}\right)}^2}$$

(6)

$$\alpha 3=\frac{{10}^{·pKa1}·{10}^{·pKa2}·{10}^{·pKa3}·\alpha 0}{{\left({10}^{·pH}\right)}^3}$$

(7)

Here α0, α1, α2, and α3 denote the relative fraction of H₃A, H₂A^-, HA^2- and A^3-, respectively. For Swedish sites, the following constants were calibrated: pK_a1 = 3.04, pK_a2 = 4.51, and pKa3 = 6.46, while the constants for Czech sites were pK_a1 = 2.5, pK_a2 = 4.42, and pK_a3 = 6.7. The median difference in fitted SD values using the two different pK_a sets is negligible, e.g., only a median of 0.88 µEq/mg C (2.7%) in the ICP Waters data set (i.e., introduced below). For consistency, we have therefore in this study only used the pK_a values for the Swedish sites.

The OAN^- is then calculated as given in Eqn. 8.

$${OAN}^{·}=TOC·\left(\frac{\propto 1·SD}{3}+\frac{2·\alpha 2·SD}{3}+\frac{3·\alpha 3·SD}{3}\right)$$

(8)

Other useful models for estimating the charge on organic matter, as documented in the scientific literature, include those proposed by Driscoll et al. [16], the ALCHEMY di- and triprotic model by Schecher and Driscoll [17], the DOM fraction model by Kortelainen [18], and the chemical equilibrium models WHAM and NICA-Donnan by Tipping [19,20] and Kinniburgh et al. [21], respectively.

3. Materials and Methods

3.1. Data Mining

The data utilized in this study are derived from several extensive monitoring programs, encompassing the ICP Waters data [22], along with the Norwegian ØkoFersk Trend [23], Reference Streams [24], and the Norwegian River monitoring program [25], as well as the 1000 Lakes study [4], commissioned by the Norwegian Environment Agency. The ICP Waters monitoring dataset, available on GitHub [26], consists of 73 656 samples collected between 1990 and 2020 from 430 sites across Europe and Northeastern America. The ØkoFersk Trend program has been conducting annual water chemistry monitoring at 78 acid sensitive Trend lakes since 1986 and weekly water chemistry monitoring, including alkalinity, at six research sites since 2001. Data from the Reference Streams monitoring program includes observations from 80 streams with different land uses sampled monthly and bi-annually since 2017. The Norwegian River monitoring program, including alkalinity since 2021, monitors monthly 22 rivers along the Norwegian coast for physical and chemical parameters. The 1000 Lakes study, conducted in the fall of 2019, sampled 2.5% of all lakes in Norway with a surface area larger than 0.04 km². The sampled lakes were distributed in a ratio of 3:2:1 among lakes from southern, middle, and northern Norway, respectively. All data for the Norwegian sites are available from the Norwegian Environment Agency [27].

3.2. Calculation of Organic Charge

The concentration of Org.^- is calculated based on charge balance principles. Below a pH of 5.5, it is assumed that there is a negligible carbonate alkalinity. However, above this threshold, HCO₃^- significantly contributes to alkalinity (i.e., α1 > 0.15). Total Alkalinity data are then utilized to estimate the concentration of weak acids, primarily HCO₃^- and Org.^-. Still, distinguishing between bicarbonate alkalinity and organic alkalinity in alkalinity measurements is challenging. In samples with pH above 5.5, the HCO₃^- concentration is calculated using Eqn. 9, derived from [28] referencing model ALK02 in [29].

$$\left[{HC{O}_3}^{·}\right]=Alkalinity·\left({{OAN}^{·}}_{SamplepH}·{O{AN}^{·}}_{pH4.5}\right)·31.62+0.2·\sqrt{Alkalinity·31.6}$$

(9)

In this context, “Alkalinity” refers to the measured total alkalinity (EN ISO 9963-1:1995), “OAN^-_{Sample pH}” represents the modelled organic charge using either the Oliver or Hruška models at the sample pH, and “ OAN^-_pH4.5” indicates the respective modelled organic charge at pH 4.5. The value 31.62 signifies the concentration of H⁺ ions at the titration endpoint (pH 4.5). All concentrations are expressed in µEq/L. The contribution of organic alkalinity to total alkalinity measurements is usually minor, encompassing solely the protonation of organic charge from sample pH to pH 4.5 [28]. Nevertheless, determining the OAN^- at sample pH and pH 4.5 to calculate the difference and thus its contribution to measured alkalinity entails a somewhat circular argument in our deductions. Additionally, it is necessary to subtract the contribution to alkalinity caused by the addition of acid to adjust the sample pH from the equivalence pH to pH 4.5.

When aiming to determine an optimal SD that establishes a 1:1 relationship between modelled OAN^- and calculated Org.^-, conceived as true value, the circular argument mentioned above poses a significant challenge. This is because increasing the SD concurrently elevates both the OAN^- and Org.^- in scenarios with carbonate alkalinity. To circumvent this issue in calculations using Eqn. 9 we fixed the SD values in the Oliver (Eqn. 3) and Hruška models (Eqn. 8) at their literature values 10.5 and 10.2 µEq/mg C, respectively. This approach ensures that the circular argument is mitigated when calculating the difference between ${{OAN}^{·}}_{Sample pH} and {O{AN}^{·}}_{pH 4.5}$ in Equation 9, thereby accounting for the contribution of organic charge to measured Alkalinity.

By estimating a concentration of bicarbonate using Eqn. 9, the Org.^- concentration can be deduced based on the anion deficit in the charge balance. If there is no anion deficit the concentration of Org.^- is assumed to be negligible and the value for Org.^- is set to zero and omitted in the regressions.

The sUVa and sVISa are the absorbency at λ254 and λ410 nm, respectively, normalized to the TOC concentration. In the Reference stream data set the sVISa is based on the Color number (i.e., Platinum number) linearly converted to absorbency at λ410 nm. The quantity of aluminum (Al) complexed to the DNOM is indirectly assessed by the concentration ratio of reactive (RAl) or directly by the ratio of non-labile Al (ILAl) to the TOC.

All calculations were conducted using Microsoft Excel, while statistical analysis was performed using MINITAB (2021).

3.3. Determining Site Density through Model Optimization

The models are calibrated by fitting the SD to achieve an optimum accuracy (i.e., OAN^- = Org.^-) by finding the SD that gives a 1:1 linear correlation anchored through origo (i.e., slope of 1) between OAN^- and Org.^-. Forcing the linear slope through origo minimizes the effect of poor model accuracy at low Org.^-.

3.4. Limitations in the Conceptual Approach

Calculation of Org.^- entails considerable uncertainties, given that it represents a small difference between the substantial sums of major cations and anions, compounded by inherent propagation of analytical errors in their constituents and in Alkalinity measurements. Furthermore, uncertainties also arise from assumptions employed to calculate HCO₃^- from alkalinity (Eqn. 9) and in the approach for assessing the protonation of organic charge from sample pH to pH 4.5. Consequently, calculated HCO₃^- concentrations can be somewhat overestimated in samples with low DNOM concentrations and may hence lead to a calculated apparent cation deficiency. In this study the concentration of Org.^- is then omitted from the regression assessment, as their concentrations are assumed insignificant and irrelevant. In assessments of median Org.^- and CD values the Org.^- concentrations were set to zero in samples with a positive charge balance.

In addition to uncertainties in Org.^- calculation, there are inherent limitations to the modelled OAN^-, as the model solely relies on a pH-governed mole fraction (α) of the CD to SD ratio, without accounting for the effects of Al and Fe complexation. Furthermore, analytical uncertainties in determining TOC and pH add to the overall uncertainty. Generally, the models are found to generate too high OAN^- at calculated Org.^- values less than 20 µEq/L.

4. Results and Discussion

4.1. Comparison of OAN^- Based on Oliver and Hruška Models

The ICP Waters dataset [22] serves as an invaluable resource for evaluating model performance across a broad spectrum of parameter values, as well as temporal and spatial dimensions. The model-derived organic charge concentrations (OAN^-) based on the Oliver (Eqn. 3) and Hruška models (Eqn. 8) are strongly correlated (R² = 0.9948, Figure 1). When employing a SD of 10.5 and 10.2 µEq/mg C, as proposed by Oliver et al. [12] and Hruška et al. [13], respectively, the regression slope is 0.80. This indicates that although the precision among these models is excellent, the accuracy is poor with the OAN^- simulated by the Oliver model being approximately 20% higher than that of the Hruška model. Decreasing the Oliver model’s SD to 8.43 µEq/mg C yields a 1:1 relationship with the Hruška model set at SD = 10.2 µEq/mg C (Figure 1). From this we infer that these models demonstrate strong compatibility, despite notable disparities in the SD values utilized to calibrate them. In this study the estimated SD for both models are therefore assessed.

4.2. Comparison of Modelled OAN^- with Calculated Org.^-

For the 1000 lakes samples [4] with total organic carbon (TOC) concentrations above 1 mg C/L the optimized SD for the Oliver and Hruška models, giving a 1:1 relationship between OAN^- and Org.^-, (Figure 2), were 19.0 and 29.0 µEq/mg C, respectively. However, these overall SD values, assessed across sites with varying SDs, appear notably high compared to the calculated average charge density (CD) values of 11.2 and 11.9 µEq/mg C. This discrepancy arises from the inclusion of many sites where the models fail to accurately simulate the data, indicated by a coefficient of determination (R²) below 0.75. Moreover, 63% of the samples had no or insignificant organic charge and thus Org.^- values were set to zero. Nevertheless, despite these discrepancies, it’s noteworthy that 71% of the Org.^- in these 1000 lakes remains unprotonated at pH 4.5. Consequently, this discrepancy does not significantly impede the comparison between the modelled organic acid anion (OAN^-) and calculated Org.^- based on charge balance, which includes HCO₃^- calculated from alkalinity using Eqn. 9.

Although the Oliver and Hruška models, when optimized for the 1000 Lakes data, yield reasonably acceptable simulations of OAN^- concentrations, explaining 65% and 72% of the variation in Org.^-, respectively, it is evident from Figure 2 that there exist subsets of sites with markedly different SDs, leading to varying OAN^- to Org.^- ratios. Clusters exhibiting high ratios likely correspond to sites where the models fail to perform effectively. To establish a more accurate estimate for SDs to be used in the Oliver and Hruška models, it becomes vital to optimize the SD based on temporal variations in TOC concentrations and pH relative to Org.^- at individual sites.

4.3. Governing Factors for Site- and Charge -Densities

Site optimized SD in the models using data from 78 Trend Lakes and 80 Reference Streams reveals an exponential increase in SD with declining TOC (Figure 3). This trend is partly mechanistically sound, as the proportion of humic to fulvic fraction generally decreases with decreasing TOC, leading to an increase in the SD of the DNOM. However, sites with a coefficient of determination (R²) dropping below 0.75 for the correlation between modelled OAN^- and calculated Org.^- exhibit large scatter in the SD with some unrealistically high SDs (e.g., > 20 µEq/mg C), particularly evident at sites with low TOC concentrations (Figure 3). Furthermore, the relationship between hydrogen ion concentration (H⁺) and SD demonstrates that at pH levels above 5.5 (i.e., H⁺ < 3.2 µEq/L), the SD in the model varies significantly. At these pH levels the levels of HCO₃^- constitutes a substantial contribution to the charge balance and, consequently, on the expense of Org.^-. Thus, the steep increase in SD at low TOC concentrations likely reflects limitations in the adopted conceptual approach. Total reactive aluminum concentrations (RAl) follow similar trends with SD as H⁺ (Figure 3), partly due to the covariance between H⁺ and RAl resulting from the pH dependence on aluminum solubility, as well as reflecting a dependence of organic-Al complexation. These anomalies lead to that the models generally tend to overestimate the OAN^- at low concentrations of Org.^- (< 20 µEq/L).

In summary, sites with low TOC concentrations inherently exhibit low OAN^-. Coupled with relatively high uncertainty in calculating low Org.^- concentrations, this results in highly uncertain optimization of Site Density (SD). Consequently, in the subsequent optimization of SDs, we have excluded sites in the linear correlations where the median TOC concentration is below 1 mg C/L and/or where the R² fails to exceed 0.75.

4.4. Spatiotemporal Variations in Functional Site Density

To establish site-specific SDs the values were optimized within the Oliver and Hruška models by fitting OAN^- to Org.^- using long-term variations in data for explanatory variables, including Alkalinity, sourced from the 78 Trend Lakes and the 6 Research Monitoring Sites in the ØkoFersk Trend monitoring program (data from [23]), as well as the 80 Reference Streams [24] and 22 major rivers draining into coastal waters [25].

After excluding sites with median TOC concentrations less than 1 mg C/L and coefficient of determination (R²) values below 0.75 in the ØkoFersk Trend lakes and Reference Stream datasets, 44 and 35 sites were retained, respectively. In Trend Lake sites where the models failed to achieve an explanatory value greater than 75%, the acid neutralizing capacity (ANC) was low or negative, indicating low or negligible Org.^- content. Conversely, Reference Stream sites where the models failed were characterized by significantly higher Alkalinity values and base cation concentrations, indicating water with high levels of bicarbonate and hard water. The median SD achieved during model optimization was 11.1 µEq/mg C for the Oliver model and 13.9 µEq/mg C for the Hruška model in the acid-sensitive Trend lakes, and 14.4 µEq/mg C for the Oliver model and 16.3 µEq/mg C for the Hruška model in Reference Streams. Over time, these SD values have however shown an increase during the monitoring period for the Trend lakes, as further assessed below. Hence, these SD values represent median spatiotemporal values.

Spatial variations in SD in the cleaned long-term dataset from Trend Lake and Reference Streams sites are primarily governed by differences in site average TOC, H⁺, and reactive aluminum (RAl) concentrations ( Figure 4 ). Data generated by the Hruška model yields similar results (not displayed). The observed increase in SD with decreasing TOC and H⁺ concentration is likely attributed to the relatively higher fraction of fulvic compared to humic moieties in the DNOM at sites with lower TOC concentrations, as elaborated below. The increase in apparent SD with decreasing RAl may partly be explained by reduced organic complexation with Al.

The spatial disparity in Oliver SD is predominantly elucidated by variations in site-specific factors such as average H⁺, RAl, and TOC concentration which reflect differences in the quality of DNOM. A multiple linear regression model incorporating H⁺ and TOC accounts for 74% and 41% of the SD variance across the 44 Trend lakes (Eqn. 10) and 35 Reference Streams (Eqn. 11), respectively.

$$SD=14.63·0.3699·H^{+}·0.2153·TOC$$

(10)

$$SD=19.20·1.723·H^{+}·0.6590·TOC$$

(11)

Here the H⁺ concentration is denoted in µEq/L, while TOC is measured in mg C/L. The p-value associated with TOC in Eqn. 10 exceeds 0.05, indicating its poor significance. This equation fails to effectively model SD beyond the threshold of 12 µEq/mg C, a limitation evident from Figs. 3 and 4. By substituting TOC with the less commonly measured RAl in the multiple linear regression analysis for the Trend lakes, a stronger equation (Eqn. 12) emerges, explaining a slightly higher proportion (78%) of the variance. Notably, all parameters exhibit significance, with p-values < 0.000. However, when applying this substitution to the data from Reference Streams, which encompass less acid-sensitive sites with a median pH of 6.4, the model’s performance does not improve.

$$SD=14.805·0.2935·H^{+}·0.01628·RAl$$

(12)

The spatial variability in SD, as derived from data collected at the 6 Research Monitoring Sites, which contain Alkalinity data since 2001, is detailed in Table 1. Among these sites, the SD was mainly governed by pH and TOC, rather than RAl. Notably, there exists a robust positive correlation between SD and the average pH of the sites (R² = 0.8295), alongside a negative correlation with average TOC concentrations (R² = 0.5964). Conversely, sulphate (SO₄^2-) concentrations demonstrate negligible explanatory power. This can be attributed to the limited number of research monitoring sites, which encompass a broad range of acid sensitivities. For instance, the Dalelva site, despite exhibiting high sulphate levels due to acid oxide pollution from nearby metal smelters in Nikel, Russia, is well-buffered, displaying high pH values and low RAl concentrations. Still, Kårvatn, having the highest SD (19.1 µEq/mg C) is a less acid-sensitive site and is exposed to the lowest levels of acid rain, reflected by high pH and low RAl, as well as low TOC (Table 1). This is in line with the fact that acid rain along with DNOM levels are important drivers for the SD.

Optimizing the SD for the modelled OAN^- to Org.^- data obtained from each of the Research Monitoring Stations yields a range of values spanning from 5 to 19 µEq/mg C (Table 1). Inherently, the Oliver model tends to exhibit limited explanatory power at Research Monitoring Sites characterized by low Org.^- levels, such as Øygardsbekken and Kårvatn (Table 1). Conversely, stronger explanatory values (R²) are observed at dystrophic sites boasting high TOC concentrations, such as Langtjern and Storgama. However, at Birkenes, despite relatively high TOC levels, the correlation remains relatively poor due to the low CD of the organic matter (2.5 µEq/mg C).

4.5. Spatial Differences in Contemporary Site Density

The SD based on recent seasonal fluctuation in monthly samples collected in 2022 from 37 Reference Streams [24] are devoid long-term temporal trends in the SD. Twenty-one streams were excluded from this analysis due to TOC concentrations less than 1 mg C/L and/or R² poorer than 0.75. Within this dataset, the median SD for the remaining 16 Reference streams, when employing the Oliver and Hruška models, stands at 13.8 µEq/mg C (STD: 4.5 µEq/mg C) and 15.8 µEq/mg C (STD: 4.8 µEq/mg C), respectively. The relatively high SD observed in this 2022 dataset, compared to the time-averaged SD from long-term monitoring data collected from the Trend lakes, aligns with a trend of increasing CD over time as assessed below. Spatial differences in SD are significantly correlated with the site-average concentrations of TOC and RAl, explaining 21% and 17% of the variation in SD, respectively (Figure 5). The explanatory power of H⁺ for spatial variation in SD is limited due to the higher pH (ranging in median values from 5.2 to 7.3) observed among these Reference Streams, which exhibit varying land use, compared to the more consistently acid-sensitive Trend lakes characterized by a median pH ranging from 4.7 to 6.6.

Among the 22 high-order rivers in the River Monitoring Program [25], coefficients of determination (R²) between OAN^- and Org.^- exceeding 0.75 were only achieved in 10 sites. Rivers that failed to fit the model generally exhibited higher pH levels (median 7.36) and lower TOC concentrations (median 2.1 mg C/L) compared to rivers that demonstrated good explanatory values (median pH 6.84 and 3.2 mg C/L). The median SD values for the Oliver and Hruška models, fitted to monthly data from 2021 to 2023, were 11.0 and 12.6 µEq/mg C, respectively. These median contemporary values are lower than those observed in the upstream Reference Streams (13.8 and 15.8 µEq/mg C for Oliver and Hruška models, respectively). Notably, as discussed below, the main difference lies in the CD, which exhibits considerably lower values in the high-order rivers compared to lower order streams, despite higher pH levels.

Samples with specific color of DNOM (sVISa) less than 0.003 are typically associated with TOC concentrations of less than 5 mg C/L [26]. Consequently, the higher SD observed at sites with low TOC concentrations (Figure 5) can be conceptually explained by the relatively lower proportion of humic compared to fulvic fractions within the DNOM. The humic fraction, characterized by a lower SD, is more susceptible to changes in solubility compared to the hydrophilic fulvic fraction as explained above. Consequently, the level of this fraction tends to be lower at sites exhibiting high ionic strength, low pH, and elevated RAl levels, commonly observed in acid sensitive regions impacted by acid rain. This conceptual framework is supported by the observation that the optimized SD decreases with increasing sVISa, indicative of a greater relative proportion of humic compared to fulvic moieties in the DNOM, within the Reference Streams dataset (Figure 6). On the other hand, a spatial Nordic study by Vogt et al. [32] noted that the specific UV absorbency (sUVa) tends to be lower in regions with heavy acid rain loading compared to less affected regions. This observation conflicts with our findings of lower medium SD at sites impacted by acid rain (conf. H⁺ and RAl vs. SD in Figure 4, and Eqns. 10 and 12, and SO₄^2- relative to SD in Table 1,) as a low sUVa implies DNOM with a relatively high fulvic acid fraction, which in turn entails a higher SD. This difference in spatial and temporal trends may partly be explained by regional differences in DNOM quality relative to temporal changes in organic Al complexation, as discussed below.

4.6. Temporal Trends in Charge Density

Conceptually, the decline in acid rain loading has an antagonistic effect on the CD, as it leads to a decrease in SD due to an increase in humic moieties, while concurrent increase in pH raises the relative protolysis (α) of the weak acid functional sites. Consequently, we would anticipate decreasing SD, while CD increases over time as acid rain diminishes. However, the decline in organic complexation of Al, an unaccounted factor in the OAN^- models, serves as a confounding factor in assessing the impact of decreasing acid rain. The significance of this is highlighted by modelling the Al complexation to fulvic acids in the Birkenes stream water using the WHAM Model VII [33]. This model gave a linear decline in organically complexed Al (NU) of 7.2 µmol Al/mg C/yr.

Simple linear temporal regression trends in CD of the DNOM at the 78 Trend Lake sites reveal a median yearly increase (slope) of 0.23 µEq/mg C/yr. Over the 34-year monitoring period, this translates to a cumulative rise of 7.81 µEq/mg C, signifying a doubling of the CD over time. In contrast, the H⁺ decreased by a median of 5.26 µEq/L from a median of 5.62 µEq/L in 1982 (equivalent to a change in pH from 5.25 in 1986 to 6.45 in 2020) in Norwegian Trend lakes. Based on the Oliver model this decrease of 5.26 µEq H⁺/L results in a median CD increase of only 0.761 µEq/mg C. Hence, the observed increase in CD cannot be attributed to enhanced deprotonation of weak acid functional groups. Since an increase in the humic to fulvic moieties of DNOM would mechanistically lead to a decrease in SD, the observed temporal rise in CD must, as also indicated above by the WHAM model, instead be attributed to reduced Al complexation with DNOM not accounted for in the model.

Spatial variations in the temporal CD slope, including samples with Org.^- set to 0, exhibit an exponential increase with decreasing ANC and TOC (Figure 7).

The pronounced temporal increases in the CD slope, exceeding 2 µEq/mg C/yr, observed at Trend Lake sites characterized by low TOC concentrations (mainly < 1 mg C/L) and negative ANC (Figure 7), may as discussed above for SD partly stem from limitations inherent in the conceptual approach. Nonetheless, uncertainties in Org.^- calculations remains low in these acidic waters, given the negligible levels of bicarbonate. The notably higher increasing slope in CD observed in lakes with negative average ANC can be largely be mechanistically explained by a temporal decrease in organic complexation with Al, as these waters with negative ANC tend to exhibit high LAl levels. Additionally, some deprotonation of weak acid functional groups with increasing pH could also contribute to these high slopes. Still, the increases in CD at these acid-sensitive sites are thus likely more substantial than those observed at more well-buffered sites. It is noteworthy that this reduced organic Al complexation, which is not accounted for in the models, implies that also the apparent SD would likewise increase significantly over time.

The intensively monitored Research Stations provide insights into both seasonal fluctuations and long-term trends at sites exhibiting variable acid sensitivity. A detailed summary of the simple temporal linear regression slopes in CD and key explanatory factors are presented in Table 2.

Changes in CD over two decades (Table 2) are at these Research monitoring sites significant relative to their average levels given in Table 1. Particularly noteworthy is the increasing trend in CD observed at the most acid rain-affected acid-sensitive sites, with the a pronounced increases noted at locations such as Birkenes, Øygardsbekken, and Storgama, which have experienced the a substantial declines in H⁺ concentrations. Interestingly, these sites also demonstrate the largest decreases in Al complexation density on the DNOM, indicated by declining ILAl/TOC ratios. Moreover, a strong negative correlation exists between the trends in CD and H⁺, explaining 58% of the variation, and the Al complexation density, which accounts for 71% of the variation in CD among the 6 Research Monitoring Stations. This increase in CD is therefore likely attributable to heightened deprotonation and reduced Al complexation resulting from the concurrent decrease in H⁺ concentrations and decline in ILAl/TOC, driven by the reduction in strong mineral acid loading (i.e., acid rain). Although the models account for the effect of changes in pH, they do not incorporate the impact of decreased Al complexation. This diminished complexation’s effect on the charged functional sites may contribute to some of the disparities between modelled OAN^- and calculated Org.^- values.

On the other hand, there has been a general concurrent temporal increase in the color of the DNOM [33], indicating a rise in the humic fraction relative to the fulvic acid fraction. This trend is likely constraining the increase in CD. Changes in climate may have had a similar effect by enhancing the influx of humic rich allochthonous DNOM to surface waters. Terrestrial material typically exhibits a higher humic character than autochthonous material [15]. This shift is likely the cause for the strong decline in CD observed at the less acid-sensitive and acid rain loaded site of Kårvatn in mid-Norway. Dalelva, a non-acid-sensitive site, received heavy acid rain loading from the now decommissioned metal smelter at Nikel in Russia. The temporal decline in CD at this site (Table 2) may be attributed to uncertainties in calculations due to low TOC concentrations and high pH, allowing for elevated bicarbonate concentrations. This is in line with the fact that the CD at both sites are positively correlated to ANC, explaining 29% of the variation.

In summary, the CD of DNOM has significantly increased over the past 20 to 40 years, especially in acid-sensitive sites experiencing a reduction in acid deposition. Moreover, the primary governing parameters for the slope in CD differ between sites. In the acid-sensitive Trend lakes, the main explanatory factors were ANC and TOC, whereas at the Research Monitoring Stations, which include both acid-sensitive and less acid-sensitive sites, the trends in CD were primarily explained by changes in acid rain loading, site acid sensitivity, and the density of Al complexed to the DNOM.

These findings have bearings for the expected future development of CD. The sUVa has been found to increase not only with decreasing acid deposition but also due to climate change [33]. This suggests that in areas where the influence of acid rain has been minimal, the SD and thus the CD has decreased with increasing DNOM, as observed in Kårvatn with low median sulphate levels (Table 1). Additionally, runoff from catchments with high biomass (forests) tends to have higher sUVa [21,31], resulting in lower SD. Furthermore, there is an effect of increasing amount and intensity of runoff, as there is more runoff through the forest floor containing DNOM with higher sUVa [31]. Consequently, it is reasonable to anticipate that the SD and thus the CD may decrease in the future due to the predicted climate change.

5. Conclusions

The Oliver and Hruška organic charge models exhibit strong explanatory value with calculated organic charge (Org.^-) based on anion deficit. Still, it is noteworthy that SD differs significantly in space and time (Table 3). This is due to disparities in the quality of the DNOM arising from differences in the relative fractions of humic and fulvic moieties. Furthermore, the charge neutralization of the anion base of these weak acid functional groups is influenced not only by pH but also by the density of aluminum (Al) complexed to the DNOM. Although the impact of pH is relatively small and is predominantly accounted for by the models, the effect of decreased metal complexation is not addressed in these models. The influence of Al complexation on the charge density (CD) of the DNOM may thus explain some of the spatial differences in SD, as well as the temporal increasing trends observed in acid-sensitive regions experiencing a decline in acid rain (Table 2).

This study highlights the complementary nature and practical utility of the conceptually grounded and empirically calibrated Oliver and Hruška models for estimating the charge concentration (OAN^-) attributable to organic anions in dystrophic freshwaters. However, calibrating and testing these models encounter limitations in sites where the concentration of organic charge (Org.^-) is low (i.e., < 20 µEq/L), such as in waters with low TOC (< 1 mg C/L), negative or low ANC, and in hard waters enriched with base cations and bicarbonate. In cases where Org.^- significantly contributes to the IB, both the site density (SD), used for model fitting, and the organic charge density (CD) values exhibit temporal changes and spatial variations.