The hydrochemical signature of incongruent weathering 1 in Iceland

11 Basaltic watersheds such as those found in Iceland are thought to be important sites of 12 CO2 sequestration via silicate weathering. However, determining the magnitude of CO2 13 uptake depends on accurately interpreting river chemistry. Here, we compile geochem14 ical data from Iceland and use them to constrain weathering processes. Specifically, we 15 use a newly developed inverse model to quantify solute supply from rain and hydrother16 mal fluids as well as allow for different mineral phases within basalts to react at differ17 ent rates, solutes to be removed via clay formation, and some Ca to be sourced from car18 bonate dissolution. While some of these processes have been considered previously, they 19 have not been considered together allowing us to newly determine their relative contri20 butions. 21 We find that weathering in Iceland is incongruent in two ways. Firstly, solute re22 lease from primary silicates is characterized by a higher proportion of Na than would be 23 expected from bulk basalts, which may reflect preferential weathering or some contri24 bution from rhyolites. This Na enrichment is further enhanced by preferential Mg and 25 K uptake by clays. No samples in our dataset (n=537) require carbonate dissolution even 26 if isotopic data (δMg, δSi, δCa, and/or Sr/Sr) are included. While some car27 bonate weathering is allowable, silicate weathering likely dominates. The complexity we 28 observe in Iceland underscores the need for inverse models to account for a wide range 29 of processes and end-members. Given that riverine fluxes from Iceland are more Na-rich 30 than expected for congruent basalt weathering, the characteristic timescale of CO2 draw31 down is likely affected. 32 Plain Language Summary 33 The chemical composition of rivers can be used to determine how processes occur34 ring at Earth’s surface affect the climate. Here, we use measurements from rivers in Ice35 land, which is an ideal test case as geologic factors are expected to simplify the inter36 pretation of river chemistry. Using a new approach, we identify a range of processes that 37 affect Icelandic rivers and quantify their relative importance. While it was assumed that 38 the results would be ”simple”, we instead find complicated behavior, which is important 39 to understand when applying similar to approaches to much larger rivers. We also find 40 that the way rivers in Iceland affect the climate is slightly different than originally ex41 pected. Nevertheless, these rivers still act to remove carbon dioxide from the atmosphere 42 and thus cool the climate. 43

in Iceland were identified through Web of Knowledge and transcribed into an Excel®database 163 through a combination of electronic transfer and manual tabulation. For this work, we 164 emphasized data from rivers, but our compilation includes a range of other water sam-165 ple types. Specifically, we distinguish ten water types: rain waters, cold springs, glacial 166 rivers, groundwater wells, hydrothermal fluids, lake waters, multi-sourced rivers, non-167 glacial rivers, soil pore waters, and unclassified rivers (Supplementary Table S1). In ad-168 dition to these classifications, we also group all of the river water data types together 169 (labeled river waters) as well as group the groundwater well, cold spring, and soil pore 170 water samples together (labeled sub-surface waters; Figure 1). We focused on collect-171 ing data for major solute concentrations (i.e., Cl, SO 4 , Na, K, Ca, Mg, Si), but if addi-172 tional trace elements and isotopic ratio data were available, we included them in our database.  When geographical coordinates of sampling locations were provided, they were con-178 verted to datum ISN1993 and recorded in decimal degrees. If a sample location was not 179 provided but a map was, an attempt was made to approximate the sample location us-  Arnórsson (1993) and Stefánsson et 183 al. (2001) are not available due to a lack of sample location map or provided locality data. 184 A map showing all of the sample locations is provided in Figure 1. 185 An effort was made to classify river water types (i.e., glacial or non-glacial river) 186 if a description was not provided but a sample location was by using visual observations 187 of river turbidity in Google Earth imagery. We note that the data reported by Gíslason 188 et al. (1996) and Moulton et al. (2000) represent average concentrations from multiple 189 measurements and that the individual data points were not provided by the authors. In 190 some cases, time-series of multiple observations from an individual site were reported (e.g., 191 Eiriksdottir et al., 2013) and, in our analysis, each sample is considered independently. 192 Other sources report single observations from a given locality (e.g., Jacobson et al., 2015). 193 In total, 22 publications were compiled amounting to 1432 observations (Louvat et al.,     The auto-correlation inherent in normalizing by either Σ + or Σ ± as well as the unit sum constraint can impact calculations such as linear regression and PCA (Pearson, 1897; Aitchison, 1983). We note that similar issues could arise in the case where solute concentrations were normalized by a single cation (e.g., X to Na ratios). To remove this autocorrelation and enable PCA, we normalize solute concentration ratios using the centered log-ratio transformation (clr), which, here, is calculated for normalization with Σ + as: where X i /Σ ± is the ratio of element i to the sum of major cations (Na, K, Ca, and Mg) 349 and sulfate measured in a river sample (subscript riv) or assumed for an end-member 350 (subscript j) and f j is the fractional contribution of end-member j to the net budget of 351 the normalization variable. To force primary minerals to dissolve and secondary min-352 erals to precipitate, we constrain the model to only find solutions where the f j values 353 are either positive (dissolution) or negative (precipitation) depending on the end-member.

354
To add strontium isotopic ratios to our mass balance framework, we utilize the equation: 87 Sr 86 Sr riv · Sr In contrast to 87 Sr/ 86 Sr, we do expect clay formation to naturally fractionate Mg, Si, and Ca isotopic ratios. Accordingly, we incorporate these isotopic systems using a different approach. To model the isotopic fractionation associated with clay formation, we adapt the model of Bouchez et al. (2013), which relates the isotopic composition of river water to the extent of element uptake using the expression: where δ X riv is the isotopic ratio of element X measured in a river sample, δ X source is the 365 contribution-weighted average isotopic ratio of element X supplied from all solute sources, 366 F X clay is the fraction of all of element X supplied to the weathering system that is taken 367 up into secondary clays, and ∆ X is the fractionation factor for element X incorporation 368 into secondary clays. We use the Bouchez  dictions based on alkalinity fluxes. We note here that R can be calculated for both the 412 gross weathering budget (before clay mineral formation; R gross ) and net weathering bud-413 get (after clay mineral formation; R net ) and that these two measures of R need not be 414 the same. For example, if clays only remove K from solution, which is a cation derived 415 from silicate weathering and not carbonate weathering, then R net will be higher than 416 R gross . This is discussed further in Section 4.3.1.

417
To characterize the degree of primary mineral incongruence, we combine the model results for the relative contributions from the multiple silicate end-members to determine the composition of the overall silicate end-member. Accordingly, the overall silicate endmember for each simulation is calculated by summing the products of the selected X/Σ ± values for each of the individual silicate end-members and their calculated gross mixing fractions (i.e., the mixing fractions re-normalized to sum to 1). Using Na as an example, the fact that only the Na-silicate end-member contributes to the Na budget and that the value of Na/Σ ± for this end-member is 1 means that the Na/Σ ± of the overall silicate end-member can be calculated from the inversion results as To assess secondary mineral congruence, we utilize the estimated fraction of the normalization variable taken up into secondary clays (f clay ). Alone, this value represents the fraction of the observed concentration of Σ ± that was taken up into secondary phases. So, we re-normalize f clay so that it instead expresses the fraction of gross solute release taken up into secondary phases (F Σ ± clay ) where: This equation assumes that the model perfectly reproduces the river observations, but 424 can be adjusted slightly for cases where there is some mismatch between the model pre-   Further, we assume that rainwater contains negligible Si and that seawater Sr to Ca ra-433 tios can be used to estimate the rainwater Sr to Σ ± ratio based on our assumed range 434 of Ca/Σ ± . For rainwater Ca, Mg, and Sr isotopic ratios, we use seawater values assum-ing that seawater-derived aerosols are the primary source of these elements in rainwa- 436 ter.

437
For the hydrothermal end-member, we take the range defined by the measurements 438 of samples in our database with temperatures over 100 • C as these are assumed to bet-439 ter represent "pure" hydrothermal fluids. We assume that hydrothermal fluids contain  of ∆ 44 Ca from -1.5 to +0.5 to account for the present uncertainty.

486
For the composition of the silicate end-member, we use a novel approach that allows for primary mineral incongruent weathering. Specifically, we break the silicate end-488 member into four separate components: Na-silicate (Na/Σ ± = 1), K-silicate (K/Σ ± = 489 1), Ca-silicate (Ca/Σ ± = 1), and Mg-silicate (Mg/Σ ± = 1). In effect, this allows the sil-490 icate end-member to take on any composition and uses the data to constrain the best-  In the original publications from which our compilation is based, many of the gen-518 eral characteristics of the solute chemistry of natural waters in Iceland were discussed 519 including proposed differences between water types (e.g., glacial versus non-glacial rivers).

520
For completeness, we re-analyze and re-state some of these observations using our larger 521 dataset. While we will primarily focus on the river water samples, various aspects of the 522 other sample types are discussed as they relate to interpreting the river water samples.   lower bounds tend to be similar with and without isotopic constraints, but, within those 692 bounds, the shape of the distribution can be substantially different (Figure 7a,d,g,j,m).

693
For F Σ + clay , including isotopic constraints tends to decrease the median values for most 694 samples (Figure 7b). In some cases, the upper bound on F Σ + clay markedly decreases when 695 isotopic constraints are included (e.g., Figure 7e,h). Accordingly, adding isotopic con-696 straints appears to help to bound the role of cation uptake into clays. As such, the high-697 est values of F Σ + clay predicted for a sample using only major elements as constraints is likely 698 an over-estimate for this parameter. (e.g., Figure 7d,f,k), but slight differences are apparent for some samples (e.g., Figure   706 7j,l).

707
While there are coherent differences in median values of R net , F ± clay , and Na/Σ +  Conservative mixing between solutes sourced from congruent basalt dissolution and 720 atmospheric deposition fails to explain a majority of the river and sub-surface water data 721 from Iceland (Figures 2, 5). This observation has been previously stated in some of the 722 original literature from which our compilation is based (e.g., Gíslason et al., 1996), al-beit sometimes using different wording and/or data analysis techniques. Nevertheless, 724 multiple competing mechanisms exist to explain this discrepancy between the null hy-  Here, we also consider the potential for individual minerals to react at distinct rates yield-733 ing a different geochemical signature than bulk basalt (primary mineral incongruence).

734
To explore the different hypotheses and help discern their relative importance, we first 735 analyze the data using mixing diagrams and then test the resulting inferences using the 736 solute mass balance model. It is important to use these two steps together as individ- ing to contribute more to the total cation budget, the relative importance of incongru-746 ent weathering must be increased. Below, we describe the evidence for each of these find-747 ings, discuss the apparent differences between water sample types (i.e., glacial versus non-748 glacial and river vs. sub-surface), and provide potential mechanisms underlying the ap-749 parently incongruent weathering taking place in Iceland.  In principle, the extent of Mg-depletion and Na-enrichment could be due to con-766 tributions from hydrothermal fluids (Figure 2a), though this would not explain the ob-767 served extent of Si depletion as most hydrothermal fluids have high Si/Σ + (Figure 2b).

768
While such a scenario appears plausible when looking at individual mixing diagrams, vi-769 sualizing the data using PCA more clearly shows that the extent of Na-enrichment in 770 the river water data is unlikely to be driven exclusively by hydrothermal contributions 771 (Figure 5a). Relative to conservative mixing between rain water and bulk basalts, the  While clay mineral formation acts to increase riverine Na/Σ + , it is likely insuffi-811 cient to explain the full extent of Na enrichment evident in the river water data such that 812 an additional mechanism for Na enrichment is required. In part, this is because the mag-  Na-rich lithologies and/or minerals could explain the general enrichment in Na/Σ + of 853 the overall silicate end-member determined from the river water data (Figure 9).

854
If the Na-enrichment of the silicate end-member is due to the preferential weath-    The results for R net from the solute mass balance model imply that the dissolu-885 tion of carbonate is not required to explain the elemental and isotopic composition of 886 Icelandic rivers (Figure 10a). Specifically, using the 5 th percentile of R net as a guide for the minimum carbonate contribution, the largest value we find for any river sample is surements, but we suggest that this would be a useful goal for future work.

929
While adding Mg, Si, Ca, and Sr isotopic constraints helps to bound cation uptake 930 into clays (Figure 7e,k), we find that C isotopic ratios help to bound possible carbon-  As an orthogonal approach to further constrain R net , we take advantage of the fact   To make a simplified comparison to previous work, we note that Gíslason et al. (1996) 1032 report that K is more mobile than Ca in Icelandic watersheds, which is not supported 1033 by our analysis (Figure 8c,d). Similarly, Gíslason et al. (1996) suggest that up to 90% 1034 of Ca is sequestered into clays, which is a much higher proportion than our upper bound   The results of the solute mass balance model are uncertain, but the few samples 1048 with the highest median carbonate contributions are from glacierized catchments (Fig-ure 13a,b). Though this is apparently consistent with the hypothesis of Jacobson et al. directly including C isotopic ratios also adds computation expense as the C isotopic ratios of weathering end-members vary independently of their X/Σ ± values thus requiring more random draws to find acceptable solutions.
For the purposes of approximating the DIC budget, we equate the positive alkalinity provided by each end-member to its DIC concentration: This assumes all DIC is speciated as HCO − 3 , which is reasonable given observations of riverine pH values. Using Monte-Carlo simulation, we take the apportionment of DIC budget between silicate, carbonate, rain water, and hydrothermal fluids for each simulation from each sample and compute a range of possible riverine C isotopic ratios given the a priori ranges in the different end-member C isotopic ratios using the equations: δ 13 C riv,pred = 4 j=1 f DIC,j · δ 13 C j (2) and 14 C riv,pred = 4 j=1 from MEANDIR f DIC,j · 14 C j random variable (3) The output of these equations can then be compared with the actual C isotopic measurement, which avoids having to require the model to exactly match the observations given that, for example, measured δ 13 C values may be elevated as a consequence of CO 2 degassing.
For the purposes of modelling C isotopic ratios, we assume that any DIC from rainwater has a δ 13 C of -7 and 14 C abundance of 105 percent modern carbon (pMC) following For DIC from hydrothermal fluids, we assume a uniform distribution of δ 13 C values between -18.8 and +4.4 following Barry, Hilton, Füri, Halldórsson, and Grönvold (2014) and a 14 C of 0 pMC.
Based on measurements of Icelandic spar from Landis (1983) and (Smalley et al., 1989), we assume a uniform distribution of carbonate δ 13 C values from -5 to -3 and a 14 C abundance of 0 pMC. For carbonate weathering driven by carbonic acid, half of the DIC will derive from the carbonic acid such that the overall values of δ 13 C and 14 C will not match pure carbonates: Here, we assume that carbonic acid derives from a mixture of atmospheric CO 2 (δ 13 C = -7 ; 14 C = 105 pMC) and CO 2 from soil respiration. Based on direct measurements of soil organic matter from Icelandic river catchments (Torres et al., 2020), we assume a uniform distribution of δ 13 C values between -32 and -24 and a uniform distribution of 14 C abundances from 60 to 130 pMC for the soil respiration end-member.
Similar to how to model carbonate weathering, we assume that DIC from silicate weathering has C isotopic compositions given by a mixture of atmospheric CO 2 and CO 2 derived from soil respiration. We draw the end-member C isotopic ratios for silicate weathering separately from carbonate weathering, which in effect allows for different CO 2 sources to drive the weathering of different mineral phases.