Estimating mineral surface area and base cation weathering rates of Spodosols under forest in British Columbia, Canada

Geoderma Regional 20 (2020) e00247

Contents lists available at ScienceDirect

Geoderma Regional journal homepage: www.elsevier.com/locate/geodrs

Estimating mineral surface area and base cation weathering rates of Spodosols under forest in British Columbia, Canada Patrick A. Levasseur a,⁎, Shaun A. Watmough b, Julian Aherne b, Colin J. Whitfield c, M. Catherine Eimers b a b c

Environmental and Life Sciences Graduate Program, Trent University, 1600 West Bank Dr., Peterborough, Ontario K9J 7B8, Canada School of the Environment, Trent University, 1600 West Bank Dr., Peterborough, Ontario K9J 7B8, Canada School of Environment and Sustainability, and Global Institute for Water Security, University of Saskatchewan, 11 Innovation Blvd, Saskatoon, Saskatchewan S7N 5C8, Canada

a r t i c l e

i n f o

Article history: Received 26 August 2019 Received in revised form 15 October 2019 Accepted 1 November 2019 Available online xxxx Keywords: Weathering rates Particle size Sulphur deposition Mineral soil Surface area

a b s t r a c t Soil mineral surface area is regarded as a key uncertainty in the estimation of base cation weathering rates, yet is rarely measured. Acidification studies rely heavily on pedotransfer functions (PTFs) that use widely available soil data to estimate mineral surface area. This study examined the relationship between soil properties and mineral surface area in soils (n = 25) from Kitimat, British Columbia, an area that is receiving elevated sulphur (S) deposition due to recent modernization of an aluminum (Al) smelter. Mineral surface area was measured on bulk soil samples using BET (Brunaeur, Emmett and Teller) gas-adsorption. Previously published particle size-based PTFs were a poor predictor of surface area in Kitimat soils (R2 between 0.42 and 0.66). Instead, mineral surface area was best predicted using a regionally-specific PTF (R2 = 0.81), which used particle size as well as the concentration of kaolinite, the most abundant clay mineral in the region. Surface area values estimated using the regionally-specific PTF were applied to the PROFILE model to calculate weathering rates for critical load estimates. These estimates predicted that none of the sites received S deposition in exceedance of their critical load for acidity. However, as surface area is largely related to kaolinite content (a mineral that does not largely contribute to weathering rates), the applicability of using surface area functions for weathering rates is questionable. Further, the texture-based PTF developed for Kitimat did not provide accurate estimates of measured surface area for other soils in Canada, particularly at surface area values exceeding 2.5 m2 g−1. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The surface area of mineral soils is influenced by several factors including particle size, organic matter (OM) content and mineralogy (Feller et al., 1992; Mayer and Rossi, 1982; White et al., 1996); consequently, mineral surface area varies widely amongst different soil types (Petersen et al., 1996). In general, soils that are sandy or contain minerals such as quartz or feldspar have a lower surface area than soils with a higher clay content dominated by 2:1 silicate clays like vermiculite (Petersen et al., 1996). Surface area affects a number of soil processes, including OM degradation (Kaiser and Guggenberger, 2000), mineral dissolution (Brantley and Mellott, 2000) and mineral weathering rates (White et al., 1996). Mineral weathering is the main long-term source of soil base cations (calcium (Ca2+), magnesium (Mg2+), potassium (K+), sodium (Na+)) in many regions and is an essential input to critical load estimates of the potential sensitivity of ecosystems to acidic deposition or timber harvesting removals (Brantley and Mellott, 2000; Hodson et al., 1996; Sverdrup and Warfvinge, 1992; Whitfield et al., 2018). ⁎ Corresponding author. E-mail address: [email protected] (P.A. Levasseur).

https://doi.org/10.1016/j.geodrs.2019.e00247 2352-0094/© 2019 Elsevier B.V. All rights reserved.

A variety of methods have been used to measure surface area, with no real consensus as to which approach best approximates field conditions (Arnepalli et al., 2008; Santamarina et al., 2002). These methods can involve wet techniques using the absorption of polar liquids and dyes (methylene blue, EGME) or dry techniques such as gas adsorption. Brunauer, Emmett and Teller (BET) gas adsorption theory is considered the foundation of surface area measurement (Brunauer et al., 1938), and many mineral weathering studies use BET gas adsorption as the primary surface area measurement tool (Hodson et al., 1998; Sverdrup and Warfvinge, 1995; Warfvinge and Sverdrup, 1992; Whitfield and Reid, 2013). The BET isotherm utilizes the pressure and volume of the adsorbent (typically N2 gas) that is needed to form a monolayer covering the entire surface of the mineral. The BET method requires OM to be removed prior to analysis as OM binds together clay particles and occludes pores thereby decreasing the amount of mineral surface available for gas adsorption (Feller et al., 1992). The BET and OM removal procedure is time-consuming and expensive and because of this, many studies have looked for alternative ways to estimate surface area using more readily available soil information (Farrar and Coleman, 1967; Hodson et al., 1998; Mayer and Rossi, 1982). Pedotransfer functions (PTFs) use more readily available soil parameters, like particle size classes to estimate surface area. One of the

2

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247

most frequently applied PTFs was developed for Swedish podzols by Sverdrup and Warfvinge (1995) and is commonly used in critical load models in Europe, Asia and North America (Guo et al., 2014; Mongeon et al., 2010; Phelan et al., 2014; Sverdrup and Warfvinge, 1992; UNECE (United Nations Economic Commission for Europe), 2004; Whitfield et al., 2018). The Sverdrup and Warfvinge model (referred to as the Sverdrup model) uses particle size classes to predict mineral surface area:
Surface Area m2 g −1 ¼ 0:08  percent clay þ 0:022  percent silt þ 0:003  percent sand

ð1Þ

Studies by Hodson et al. (1998) and Whitfield and Reid (2013) determined that this particle size-based function did not accurately estimate surface area when applied to soils from different regions (Scotland, UK and Alberta, Canada, respectively). Both studies suggested that consideration of mineralogy in addition to particle size was necessary to improve the accuracy of PTF predictions. Despite these concerns the Sverdrup model is still widely applied to predict surface area for critical load calculations across the globe. Mineral surface area and base cation weathering rates are key inputs to critical load models that predict the potential impact of acidic deposition. While sulphur (S) deposition decreased in Canada by approximately 60% between 1990 and 2012 (UNECE (United Nations Economic Commission for Europe), 2016), S deposition has increased in some regions of Canada as a result of local point sources. Specifically, in 2012, an aluminum (Al) smelter in the largely pristine forested region of Kitimat, British Columbia, Canada was permitted to increase S emissions by over 50% (27 t day−1 to 42 t day−1) (ESSA Technologies et al., 2013). Before beginning operations, the company was required to evaluate impacts on all ecosystems and populations that were expected to receive elevated S deposition (ESSA Technologies et al., 2013). The ESSA Technologies et al. (2013) study utilized mineral surface area predictions from the Sverdrup PTF as inputs to calculate critical loads for the region and found that critical loads were generally high and only soils directly around the smelter would receive S deposition in exceedance of the critical load (ESSA Technologies et al., 2013; Williston et al., 2016). However, the study raised concerns over uncertainties in mineral weathering rate estimates and the use of the Sverdrup PTF (ESSA Technologies et al., 2013). Thus, the objectives of this study were to a) establish relationships between measured mineral surface area and more widely available soil properties, b) develop a PTF for the Kitimat region based on mineral surface area measurements, c) evaluate regional PTF versus previously published PTFs on soils from Kitimat and across Canada, d) use the best-suited PTF to calculate critical loads for S deposition in Kitimat. Mineral surface area can be predicted from easily measured soil properties but a locally-developed PTF would provide a better estimate of surface area than previously published models. 2. Methods 2.1. Kitimat study area and soil collection The projected area of enhanced S deposition is approximately 2895 km2 and extends from the mouth of the Douglas Channel, just south of the smelter in Kitimat to the city of Terrace, British Columbia in the north (Fig. 1). This enhanced deposition zone includes the entirety of the Kitimat Valley as well as the surrounding mountains, and encompasses an elevation range of 300–1700 m. The ocean and mountain influences make the entire affected region one of the wettest places in Canada, with average precipitation and temperature (values from 1981 to 2010) of 1781 mm and 5.6 °C, respectively (Environment Canada, 2015). The enhanced deposition zone is in the Coastal Western Hemlock ecozone and is dominated by conifers such as western hemlock (Tsuga heterophylla), balsam fir (Abies balsamea), sitka spruce (Picea sitchensis) and red cedar (Thuja plicata). The dominant soils in

the area are Ferro-Humic Podzols (Humic Cryorthods in U.S. soil taxonomy) (Soil Classification Working Group, 1998), the cold and wet nature of these soils leads to a thick forest floor, high OM concentrations in mineral horizons and moderately acidic soils (Table 1). Twenty-five soil sampling sites were selected to encompass the range in parent materials in the region (ESSA Technologies et al., 2013), including bedrock (5 sites), morainal (5), glaciofluvial (6), marine (3), fluvial (5) and colluvium (2). At each site, a 10 m × 10 m quadrat was established, and mineral soil was sampled using a soil auger at each corner and at the center point of each plot at three depths (0–10 cm, 15–25 cm, 40–50 cm) after first removing the forest floor. Soil samples were sampled by depth versus horizon to standardize procedures between field teams. Soil cores for the determination of bulk density were taken using a bulk density hammer (the height and diameter of the core was 5 cm) at each depth at the centre of the plot. This field sampling occurred in the summer of 2012. 2.2. Laboratory analysis 2.2.1. Physical analyses All soils were oven dried at 105 °C and crushed with a mortar and pestle prior to being sieved to b2 mm, samples from each site were composited by depth, weighted by bulk density. Sites were composited by depth to determine characteristics of entire soil profile. Specifically, for each site the bulk density at each depth was divided by the average bulk density for all depths to determine the percentage of soil from that depth to be put into the composite. Loss on ignition (LOI) was performed to determine OM concentration by ashing oven-dry soils for 8 h at 475 °C in a Fisher Scientific Isotemp Muffle Furnace. Samples were stored in a desiccator and weighed before and after ignition. Ashed soil samples were then used for particle size analysis using a Horiba Particle size Analyser LA-950 V2 to measure the percent sand, silt and clay. Each sample was analyzed in duplicate, suspended in a solution of 40 g L−1 of sodium hexametaphosphate to reverse osmosis (RO) water. Laser diffraction uses a volume-based approach to determine particle size. Mineralogy was retrieved from ESSA Technologies et al. (2013). All mineralogy values reflect the percentage of the nonquartz minerals that they make up. Quartz does not play a role in mineral weathering and has been removed from mineralogy results (ESSA Technologies et al., 2013). Soil pH was measured in a 1:4 soil: deionized water slurry. Cation exchange capacity (CEC) was measured on 5 g of dried soil extracted with unbuffered ammonium chloride (NH4Cl) and sodium chloride (NaCl) (Lozano et al., 1987). Cation exchange capacity extracts were analyzed for NH4 by colourimetry with a Bran and Luebbe Autoanalyzer 3. 2.2.2. Surface area measurement Organic matter and amorphous iron (Fe) and aluminum (Al) oxides may coat and cement individual soil particles and thereby alter their reactive surface area (Feller et al., 1992; McKeague and Day, 1966), and because of this both OM and Fe and Al oxyhydroxides were removed from soils prior to BET surface area analysis. Organic matter was removed using sodium hypochlorite (NaOCl) adjusted to pH 8 using hydrochloric acid (HCl), the NaOCl method has been reported to effectively remove OM without altering the structure of soil minerals (Mikutta et al., 2005; Kaiser and Guggenberger, 2003). Falcon tubes with 50 mL of NaOCl solution and 1 g of soil were hand shaken until soil particles that had been bound together dissolved. These falcon tubes were then placed vertically on a shaker table for 12 h at 60 oscillations per minute and centrifuged at room temperature at 2000 rpm for 20 min. After centrifuging, the solution and all suspended OM were decanted from samples. Due to the organic nature of these soils, this oxidation process was repeated five times (Kaiser and Guggenberger, 2003). After the fifth NaOCl rinse the same falcon tube was filled with 50 mL of RO water and then hand shaken until bound

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247

3

Fig. 1. The location of the 25 soil sampling sites in the Kitimat region. Sites were selected to cover the range of parent materials within and outside of the zone of projected increased S deposition. The purple star represents the aluminum smelter.

soil particles dissolved and centrifuged at room temperate at 2000 rpm for 20 min. Following NaOCl treatment, soils were treated with 40 mL of a 4:3 solution of 2 M acid ammonium oxalate AAO and oxalic acid to remove Fe and Al oxyhydroxides. Tubes were hand shaken until bound soil particles dissolved and then placed on a shaker table at 60 oscillations per minute for 8 h in the dark. Samples were then centrifuged at room temperature at 2000 rpm for 20 min, decanted, filled with 50 mL of RO water, hand shaken until bound soil particles dissolved and centrifuged again at room temperature at 2000 rpm for 20 min before decanting once more. Treated soils were then oven-dried at 40 °C for 48 h before BET analysis. Table 1 Average, maximum and minimum values of soil characteristics for bulk density weighted averages of three depths (0–10 cm, 10–25 cm and 25–50 cm)at 25 sites where BET surface area was measured. Soil characteristic

Average Minimum Maximum Standard deviation

Sand (%) Silt (%) Clay (%) OM (%) Elevation (masl) Bulk Density (g cm−3) pH CEC (mmolc kg−1) BET Surface Area (m2 g−1) Apatite (%) Calcite (%) Fe-Chlorite (%) Hornblende (%) K-Feldspar (%) Kaolinite (%) Muscovite (%) Plagioclase (%)

63.9 32.8 2.8 8.6 344.9 0.8 5.1 114.2 2.2 0.7 0.2 10.7 9.8 24.5 14.0 9.8 30.3

33.0 8.8 0.0 2.1 6.0 0.4 4.6 33.6 0.1 0.1 0.0 0.7 3.0 10.9 0.3 1.1 14.1

91.7 61.6 6.4 19.3 1249.0 1.3 6.1 160.7 6.9 1.7 2.2 23.5 19.4 65.2 30.6 16.7 63.8

6.9 6.4 0.6 3.4 157.7 0.2 0.4 25.4 2.2 0.4 0.6 5.4 4.5 12.4 10.1 4.0 12.4

The BET nitrogen gas adsorption isotherm was used to measure mineral soil surface area following methods described in Brunauer et al. (1938) using a Micromeritics™ Gemini VII analyser. In brief, Micromeritics tubes were degassed with N2 using a Micromeritics ™ Flow Prep Degasser for 20 min at 30 °C and weighed on an analytical balance. The tubes were then filled with 2 g of sample and degassed for another 24 h at 30 °C before re-weighing. The dewar was filled with 500 mL of liquid nitrogen (LN2) and left to equilibrate with the temperature of the dewar, and the dewar was refilled every three samples as LN2 evaporated. Measurement of initial pressure (P0) was taken with an empty tube before sample analysis and sample tubes were filled with Micromeritics 3 mm glass beads to minimize free space. BET measurements were run on the 25 sites in triplicate.

2.3. Pedotransfer functions Relationships between BET measured surface area and site elevation, pH, bulk density, mineralogy, OM concentration, cation exchange capacity and particle size were tested for the 25 study sites through Pearson correlation. All significantly correlated variables were applied to a stepwise linear regression to develop a Kitimat specific PTF. A PTF including only particle size classes was also developed for Kitimat soils using stepwise linear regression. While PTFs using mineralogy and particle size have been found to better predict surface area in their respective regions, adding mineralogy limits the transferability of these PTFs (Whitfield and Reid, 2013). The surface area estimates from this Kitimat particle size model (referred to as the texture-based model) were compared with estimates from other previously established PTFs developed by Sverdrup and Warfvinge (1995). Hereafter, the PTF models published by Hodson et al. (1998) and Whitfield and Reid (2013) will be referred to as the Hodson and Whitfield functions, respectively.

4

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247

Chlorite were the most abundant clay minerals and averaged 14.0% and 10.7% of mineral composition respectively (Table 1).

2.4. Weathering rates and critical loads in Kitimat As nitrogen (N) deposition to the region was estimated to be very low (ESSA Technologies et al., 2013) the focus of critical load determinations was on S (CL(S)). Atmospheric deposition inputs of base cations (BCdep) and Cl were not available for the study area, and so a conservative approach was taken in assuming that base cation and Cl deposition is 0 (eq. 2). CL ðSÞ ¼ BC w –Bcu −ANC leðcrit Þ

ð2Þ

Base cation uptake (Bcu = Ca + Mg + K) values were estimated using a combination of estimated biomass removal rates in timber harvest and literature values of the concentration of base cations within that biomass (ESSA Technologies et al., 2013). Mineral weathering rates were calculated using PROFILE version 5.1 (Sverdrup and Warfvinge, 1988). PROFILE includes individual variables such as soil moisture, soil bulk density, runoff, soil temperature and mineralogy. Weathering rates (BCw) were calculated assuming a soil depth of 0.5 m. This approach was used to approximate how altering surface area would affect weathering rates for individual soil samples. Critical acid neutralizing capacity leaching (ANCle(crit)) was calculated using eq. 3. 2

ANC leðcrit Þ ¼−Q 3 ð1:5 BC w þ BC dep – Bcu =ðBc : AlÞ K gibb
–1:5 BC w þ BC dep – Bcdep =ðBc : AlÞ



13

ð3Þ

Long-term runoff (Q) was calculated following Moore et al. (2012). The ratio Bc:Al was set to 1.0, which is widely used in other studies (Cronan and Grigal, 1995). The gibbsite equilibrium constant (Kgibb) was set to a regional default of 8.5 (ESSA Technologies et al., 2013). Critical load exceedances were calculated by subtracting S deposition from CL(S). Sulphur deposition data were obtained from ESSA Technologies et al. (2013) based on CALPUF (a Langrangian puff modelling system) modelled Kitimat deposition using an average of meteorological years of 2006, 2008 and 2009. Weathering rates and critical loads were calculated for 24 sites. Deposition data were not available for all 25 study sites (Fig. 1) and so weathering rates and critical loads were calculated for 18 of the original 25 study sites plus an additional six sites where S deposition and particle size data were available. 2.5. Statistical analysis

3.2. Physicochemical properties relationship to surface area Measured surface area varied greatly amongst the 25 sites, ranging from 0.04 to 6.86 m2 g−1. The surface area of soils derived from bedrock parent material was significantly lower than almost all other surficial materials including morainal (p b .05), fluvial (p b .05), glaciofluvial (p b .01) and marine (p b .005; Fig. 2). There was no significant correlation between surface area and site elevation (r = −0.17), OM (r = −0.04), bulk density (r = −0.003), pH (r = 0.16), hornblende (r = −0.32), KFeldspar (r = −0.28) or CEC (r = 0.21). Surface area was significantly positively correlated with percent clay (r = 0.88), percent silt (r = 0.70), Fe-chlorite (r = 0.69) and kaolinite (r = 0.85) and significantly negatively correlated with percent sand (r = −0.75) and plagioclase (r = −0.66) (Fig. 3). 3.3. Comparison of texture-based PTF estimates to measured BET values A PTF using only soil particle size was created for soils from the Kitimat region (Table 2). A comprehensive PTF was also developed including all variables significantly correlated with mineral surface area. While plagioclase and Fe-chlorite were found to be significantly correlated with mineral surface area their inclusion in the comprehensive PTF increased AIC scores. The final comprehensive model included the three particle size classes (sand, silt and clay) as well as kaolinite, the most abundant clay mineral in the Kitimat region (Table 2). The texture and comprehensive models yielded similar results suggesting there was little benefit of including clay type (Table 2). In order to evaluate the transferability of previously developed texture-based PTFs to the Kitimat region, the model predictions of Sverdrup, Hodson and Whitfield were compared with BET measured SA of Kitimat soils. The Sverdrup (Adjusted R2 = 0.47) and Hodson (Adjusted R2 = 0.42) functions had very poor fits, while the Whitfield (Adjusted R2 = 0.66) PTF was marginally better. The NRMSE results showed a similar trend with error as the Sverdrup (NRMSE = 1.04) and Hodson models (NRMSE = 1.03) had much larger degrees of error compared to the Whitfield (NRMSE = 0.71) or Kitimat model (NRMSE = 0.44). Inaccuracy of surface area estimates increased with measured values for all of the Sverdup, Hodson and Whitfield models, surface area estimates underpredicted measured values N3 m2 g−1 by up to 400%.These results indicated a locally-derived PTF was necessary to accurately predict BET measured surface area (Fig. 4).

Normality of all data was tested using the Kolmogorov-Smirnov test and non-normal variables were log or square root transformed. If variables were still non-normal, non-parametric tests were used. Differences in surface area values between surficial materials were assessed using a Mann-Whitney U test, these differences were considered significant at p b .05. Pedotransfer functions were determined using a backward stepwise linear regression and were evaluated using Akaike Information Criterion (AIC) scores, linear regression, bootstrapping and normalized root-mean-square error (NRMSE). Differences in weathering rate estimates were assessed using a Mann-Whitney U test, these differences were considered significant at p b .05. All statistical analysis was done on R version 3.3.0 (Rstudio Team, 2015). 3. Results 3.1. Site characteristics Mineral soils in the Kitimat region were generally sandy (average 64%), acidic (pH b 6.2) and high in OM concentration (Table 1). Plagioclase was the most abundant mineral (non-quartz fraction) averaging 30.3% followed by K-Feldspar (average 24.5%). Kaolinite and Fe-

Fig. 2. BET measured surface area (m2 g−1) for all known surficial materials in the study region. Bars represent average surface area. Error bars are standard error amongst surficial geologies and numbers are sample size. Letters represent significant differences.

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247

5

Fig. 3. Square root BET measured surface area (m2 g−1) compared against other soil characteristics. The characteristics that were significantly correlated with square root BET were clay (%), square root plagioclase (%), Fe-Chlorite (%) and kaolinite (%). Shaded area represents 95% confidence interval.

3.4. Texture-based PTFs applied to soils from across Canada One factor that might influence the transferability of texture-based PTFs from one site to another is the method of texture analysis. In order to compare the three previously published PTFs and the newly developed Kitimat PTF, soils were collected from 19 sites across Canada and analyzed for soil texture and BET surface area using consistent methods. All four PTFs had very low R2 values (0.08–0.15) overall inaccuracies appeared to increase with larger surface area values as surface areas higher than 3 m2 g−1 were greatly underestimated, sometimes by as much as 4.5 m2 g−1 (Fig. 5).

3.5. Weathering rates, critical loads and exceedances Weathering rates calculated with PROFILE using surface area estimates from the Sverdrup, texture and comprehensive PTF were compared for all 24 sites (Table 3). Average weathering rates varied between 1.25 and 2.53 kmolc ha−1 yr−1. Weathering rates estimated using the Sverdrup PTF were significantly lower than weathering rates estimated using both the comprehensive and texture-based PTFs (p b .05). There was no significant difference in weathering estimates between the texture or comprehensive PTFs (V = 161, p = .5). Critical loads were estimated for 24 sites using surface area values from the texture PTF. These critical load estimates ranged from 1.13 to 15.65 kmolc ha−1 yr−1, with an average of 7.11 kmolc ha−1 yr−1

(Fig. 6). Only three of the 24 sites had critical loads below 2.00 kmolc ha−1 yr−1, all three sites were found on bedrock parent material. Out of the 24 sites, none received S deposition in exceedance of their critical load. 4. Discussion 4.1. Mineral surface area and other physicochemical soil properties BET-measured surface area values for the Kitimat soils ranged from 0.04 to 6.84 m2 g−1. This range is similar to BET-measured values found in other studies. Lower values (ranging from 0.54 to 2.26 m2 g−1) were found in podzols and gleysols of Scotland that were developed on bedrock of low weathering material (granite, schists and gneiss) (Hodson et al., 1998). Higher values (range of 0.1 to 14 m2 g−1) were found in brunisols and luvisols in Alberta with a higher clay content (median 5.01%) (Whitfield and Reid, 2013). Soils from areas dominated by bedrock had a significantly lower surface area (average of 0.75 m2 g−1) than all other surficial materials (morainal, fluvial, glaciofluvial and marine). Soils formed on bedrock-dominated areas that have undergone little weathering are expected to have a much lower surface area than soils derived from finer surficial materials. Previous studies have found that soil cation exchange capacity (Farrar and Coleman, 1967), soil particle size (White et al., 1996) and mineralogy (Hodson et al., 1996) are correlated with mineral surface area. Out of the 12 variables tested in this study, clay (+), silt (+),

Table 2 Equations for regionally developed PTFs and their evaluation scores (including Akaike Information Criterion (AIC), adjusted R2 (Adj R2) and normalized root mean square error (NRMSE)) for estimating surface area. Coefficients Model

Intercept

Clay (%)

Silt (%)

Sand (%)

Kaolinite (%)

AIC

Adj R2

NRMSE

Texture Comprehensive

9.29 6.72

0.30 0.20

−0.10 −0.07

−0.09 −0.06

0 0.02

−50.19 −51.95

0.80 0.81

0.44 0.47

6

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247

Fig. 4. Texture-based pedotransfer functions (PTF) estimates compared with BET measured mineral surface area (BET SA). Lines represent 1:1 ratio. Results demonstrate that previously developed PTFs (Sverdrup, Hodson and Whitfield) are not transferable to the Kitimat region. Adjusted R2 (Adj R2) and normalized root-square mean error (NRMSE) are used to assess fit of model.

sand (−), plagioclase (−), Fe-chlorite (+) and kaolinite (+) were the only soil parameters that were significantly correlated with BET surface area. The lack of relationship between CEC and surface area could be explained by the type of clay minerals that are present in these soils. The strongest relationships between CEC and surface area were found in soils that had three major clay minerals: kaolinite, illite and montmorillonite (Farrar and Coleman, 1967). Similar to the present study, Petersen et al. (1996) found that Danish soils with a wide variety of clay minerals did not have a significant relationship between surface

area and CEC. Many studies have found a strong relationship between surface area and texture, as smaller particles such as clay have larger surface areas (Mayer and Rossi, 1982; White et al., 1996). Plagioclase, a mineral with low surface area (White and Brantley, 2003), is the most abundant mineral in the Kitimat soils (excluding quartz) and is negatively correlated with surface area. Kaolinite was the most abundant clay mineral, followed by Fe-chlorite and as their abundance increased so did surface area. The relationships between plagioclase and kaolinite concentrations with surface area of bulk soil samples is

Fig. 5. Texture-based pedotranfser function (PTF) estimates compared with BET measured mineral surface area (BET SA) of soils from 19 sites across Canada. Lines represent 1:1 ratio. Adjusted R2 (Adj R2) and normalized root-square mean error (NRMSE) are used to assess fit of model.

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247 Table 3 Weathering rates estimated using PROFILE with surface area values from three pedotransfer functions (PTFs) for 24 sites where mineralogy and particle size were available. Letters in superscript denote significant differences (p b .05) between weathering rate estimates. PTF

Average (kmolc ha−1 yr−1)

Maximum (kmolc ha−1 yr−1)

Minimum (kmolc ha−1 yr1)

Sverdrup Texture Comprehensive

1.25a 2.53b 2.39b

4.70 5.50 4.86

0.17 0.31 0.31

expected as Clausen and Fabricus (2000) reported that kaolinite had a surface area of 8 m2 g−1 while White et al. (1996) measured plagioclase surface area values between 0.26 and 1.48 m2 g−1. 4.2. Pedotransfer function development The texture-based function developed by Sverdrup and Warfvinge (1995) is used to estimate surface area for calculating weathering rates and critical loads in North America and Europe (Mongeon et al., 2010; Ouimet and Duchesne, 2005; Sverdrup and Warfvinge, 1995; Whitfield et al., 2006; Williston et al., 2016). Previous studies have found that this texture function typically under-predicts measured surface area when applied to soils outside its study region (Hodson et al., 1998; Whitfield and Reid, 2013). This also proved to be true for soils from the Kitimat region. Surface area predictions from the Sverdrup function, along with the Hodson and Whitfield texture-based functions, were poorly correlated with BET measured surface area. A PTF that incorporated particle size relationships specific to the Kitimat region proved to predict surface area values that better reflected BET measurements. The explanation of why texture PTFs removed from their study region under-predict surface area largely remains unclear. This inaccuracy has previously been linked to differences in particle size and mineral

7

surface area relationships across the globe (Hodson et al., 1998; Whitfield and Reid, 2013). However, this inconsistency could also be due to the use of different techniques for particle size measurement. The current study used laser-diffraction, a volume-based approach to measure particle size. Whitfield and Reid (2013) similarly used laserdiffraction, but did not use a sodium hexametaphosphate solution to deflocculate finer particles. Hodson et al. (1998) and Sverdrup and Warfvinge (1995) used the granulometric method of sieving, a massbased approach to determine texture. Mass-based approaches to measuring texture have typically been found to have considerably higher percentages of clay fractions compared to volume-based approaches (Eshel et al., 2004). Nevertheless, when these PTFs were applied to soils from across Canada using consistent methodology to measure particle size and surface area the estimates from all four functions were not significantly correlated with BET measurements (Fig. 5). Even the Kitimat based function that was developed using the laser-diffraction particle size measurement under-predicted BET when values exceeded 2.5 m2 g−1 (Fig. 5). The under-prediction at higher surface area may be because PTFs developed in regions with predominately 1:1 clays (e.g. Kaolinite) such as in the present study do not accurately capture the higher surface areas associated with 2:1 clays (e.g. vermiculite) that may dominate in soils with higher measured surface areas. As a result, the influence of clay minerals that have large surfaces on total surface area may be lost when considering a single clay value. Incorporating mineralogy within PTFs can increase the accuracy of surface area predictions (Hodson et al., 1998; Whitfield and Reid, 2013). As different minerals weather in different patterns, adding mineralogy to a PTF provides more context on how much surface area will be exposed throughout the weathering process. Plagioclase and Kfeldspar are the main primary minerals and together make up 54.8% of the non-quartz minerals in Kitimat soils. Both minerals have a high weathering potential and contribute Na+, Ca2+ and K+ to the soil solution (Deer et al., 1992). Kaolinite, a secondary mineral sometimes formed from the weathering of plagioclase, is the most dominant clay

Fig. 6. Map of critical loads for the 24 study sites where particle size and deposition data were available. No sites received sulphur deposition in exceedance of the critical load.

8

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247

mineral in the region making up on average 14% of the non-quartz soil minerals. Kaolinite appears to be the main driver for bulk soil surface area, however kaolinite does not actively contribute base cations (Deer et al., 1992). This brings into question the appropriateness of including a generic soil surface area value that is heavily influenced by clay content for mineral weathering rates as opposed to mineralspecific surface areas. Presently however, the use of bulk surface area to calculate weathering is currently the most commonly used method, as field measurements of individual mineral surface areas are not widely available. 4.3. Critical loads and exceedances With the above limitations of measuring bulk surface area in mind, the average critical load of S at the 24 sites was 7.11 kmolc ha−1 yr−1, and none of the sites received S deposition in exceedance of the critical load regardless of which PTF was used. Critical loads estimated in this study are amongst the highest reported for Canada, including the Georgia Basin, southern British Columbia (0.94 to 4.53 kmolc ha−1 yr−1) (Mongeon et al., 2010), eastern Canada (0.56 to 2.06 kmolc ha−1 yr−1) (Ouimet et al., 2006), Manitoba (1.26 kmolc ha−1 yr−1) and Saskatchewan (1.15 kmolc ha−1 yr−1) (Aherne and Watmough, 2006). Despite significant differences between weathering rates depending on the PTF applied, critical loads of S deposition for the Kitimat region remain sufficiently high to avoid exceedance even when the lowest weathering rate estimates are applied. Similar to results from ESSA Technologies et al. (2013) and Williston et al. (2016) this study suggests that there is little risk of adverse effects of increased S deposition regardless of which surface area estimate is used. There are further additional uncertainties aside from mineral surface area estimations in these critical load calculations including: Bc:Al thresholds, timber harvesting and N deposition. These estimates use a critical Bc:Al value of 1 which is the standard for coniferous forests (Cronan and Grigal, 1995). However, some Canadian studies have found that setting a Bc:Al ratio of 1 was not effective for preventing all individual soil horizons from exceeding the toxicity threshold and recommended a Bc:Al of 10 (Ouimet et al., 2001). We also assessed the impact of increasing the critical Bc:Al to 10 and still found that none of the sites exceeded their critical loads and 9 sites had exceedances below 2.00 kmolc ha−1 yr−1 (Appendix). Timber harvesting estimates were a source of uncertainty in this study as values were a combination of base cation concentrations for tree species from literature and timber removal estimates (ESSA Technologies et al., 2013). In parts of Europe and North America, acidic deposition has been dramatically reduced in the past decade, however, an increase in base cation removal from harvesting is increasingly important with respect to soil acidification (Iwald et al., 2013; Reid and Watmough, 2016). Finally, our critical load calculations do not include N emissions for the region as current N deposition in the region is low. However, N emissions could change in the near future as liquefied natural gas (LNG) facilities have been proposed for the region potentially increasing nitrous oxides emissions from 11.9 to 30.3 t per day (Williston et al., 2016). 5. Conclusions This study found that previously published PTFs using soil particle size did not accurately predict measured surface area for soils from the Kitimat region. Applying a regional particle size class–based model predicted BET surface area measurements with a good deal of accuracy. However, when these texture–based PTFs are applied to soils across Canada, BET measured surface area is under-predicted in soils with higher surface areas (N3 m2 g−1). This may be because clay content has a strong influence on mineral surface area and not taking clay type into account may lead to inaccuracies in surface area estimates. Finally, a PTF using particle size and kaolinite was found to be the most effective model to predict surface area for the Kitimat region and suggests that

bulk surface area is influenced by a mineral that does not contribute to mineral weathering. Even though mineral surface area can be well predicted from widely available soil properties, greater consideration as to how mineral surface area values are measured and how these values are used in the estimation of weathering rates and critical loads is needed. Acknowledgements The authors would like to thank Kevin Adkinson, Phaedra Cowden, Conor Gaffney, Eric Grimm, Candace Hilchuck, Chris Jones, Liam Murray, Emily Olmstead and Carolyn Reid for their help in the field and lab. The authors would like to thank Rio Tinto BC Works, Canada for logistical and safety support during field work and financial support. The authors would also like to thank the NSERC committee for financial contributions to this project through an NSERC Discovery grant. Appendix. Critical load (CL(10)) and acid neutralizing capacity (ANC) outputs when base cation to aluminum ratio (Bc:Al) was raised to 10.

Site

CA004 CA008 G0028 GD003 GD013 GR002 GR005 LM006 LM009 OG001 OG003 QD006 QD007 QD012 QM003 S022 SSS001 SSS005 SSS006 VA001 VA006 VC002 VC003

keq ha−1 yr−1

keq ha−1 yr−1

keq ha−1 yr−1

ANC(10)

CL(10)

Exceedance(10)

−1.09 −0.97 −0.51 −0.56 −0.65 −0.77 −0.49 −0.95 −0.83 −0.72 −0.22 −0.44 −0.53 −1.23 −0.41 −0.72 −0.31 −0.27 −0.93 −0.64 −1.01 −0.40 −0.67

5.43 5.88 1.95 2.80 3.63 4.61 1.92 6.21 5.22 2.70 0.54 2.09 1.59 7.77 1.20 4.32 0.76 0.97 5.80 3.57 6.26 1.66 4.30

3.95 5.59 1.78 2.51 3.42 4.40 1.87 5.97 4.99 2.65 0.48 1.74 1.51 7.42 1.14 4.17 0.72 0.91 5.64 3.43 6.13 1.59 4.05

References Aherne, J., Watmough, S.A., 2006. Calculating critical load of acid deposition for forest soils in Manitoba and Saskatchewan. Final Report. 2006. Canadian Council of Ministers of the Environment 67pp. Arnepalli, D.N., Shanthakumar, S., Hanumantha, R., Singh, D.N., 2008. Comparison of methods for determining specific-surface area of fine-grained soils. Geotech. Geol. Eng. 26, 121–132. Brantley, S.L., Mellott, N.P., 2000. Surface area and porosity of primary silicate minerals. Am. Mineral. 85, 1767–1783. Brunauer, S., Emmett, P.H., Teller, E., 1938. Adsorption of gases in multimolecular layers. J. Am. Chem. Soc. 60, 309–319. Environment Canada, 2015. Canadian Climate Normals. p. 1981e2010. Clausen, L., Fabricus, I., 2000. BET measurements: outgassing of minerals. J. Colloid Interface Sci. 227, 7–15. Cronan, C.S., Grigal, D.F., 1995. Use of calcium/aluminum ratios as indicators of stress in forest ecosystems. J. Environ. Qual. 24, 209–226. Deer, W.A., Howie, R.A., Zussman, J., 1992. An Introduction to Rock-Forming Minerals. second ed. Longman, New York. Eshel, G., Levy, G.J., Mingelgren, U., Singer, M.J., 2004. Critical evaluation of the use of laser diffraction for particle-size distribution analysis. Soil Sci. Soc. Am. J. 68, 736–743. ESSA Technologies, Laurence, J., Limnotek, Risk Sciences International, Rio Tinto Alcan, Trent University, Trinity Consultants, University of Illinois, 2013. Sulphur Dioxide Technical Assessment Report in Support of the 2013 Application to Amend the P200001 Multimedia Permit for the Kitimat Modernization Project. 2 p. 450 Final Technical Report. Prepared for Rio Tinto Alcan, Kitimat, B.C.

P.A. Levasseur et al. / Geoderma Regional 20 (2020) e00247 Farrar, D.M., Coleman, J.D., 1967. The correlation of surface area with other properties of nineteen British clay soils. J. Soil Sci. 18, 119–124. Feller, C., Schouller, E., Thomas, F., Rouiller, J., Herbillon, A.J., 1992. N2- BET specific surface areas of some low-activity clay soils and their relationships with secondary constituents and organic matter contents. Soil Sci. 153, 293–299. Guo, P., Wang, Y., Wang, Y., Wang, R., Hu, B., Tang, X., 2014. Evaluation of critical Sulphur loads based on weathering rate modeling in three gorges reservoir area, China. Vegetos. 27, 82–89. Hodson, M.E., Langan, S.J., Wilson, M.J., 1996. A sensitivity analysis of the PROFILE model in relation to the calculation of weathering rates. Appl. Geochem. 11, 835–844. Hodson, M.E., Langan, S.J., Meriau, S., 1998. Determination of mineral surface area in relation to the calculation of weathering rates. Geoderma. 83, 35–54. Iwald, J., Löfgren, S., Stendahl, J., Karltun, E., 2013. Acidifying effect of removal of tree stumps and logging residues as compared to atmospheric deposition. For. Ecol. Manag. 290, 49–58. Kaiser, K., Guggenberger, G., 2000. The role of DOM sorption to mineral surfaces in the preservation of organic matter in soils. Org. Geochem. 31, 711–725. Kaiser, K., Guggenberger, G., 2003. Mineral surfaces and soil organic matter. Eur. J. Soil Sci. 54, 219–236. Lozano, F.C., Parton, W.J., Lau, J.K.H., Vanderstar, L., 1987. Physical and Chemical Properties of the Soils at the Southern Biogeochemical Study Site. Dorset Research Centre, Ontario Ministry of the Environment, p. 76. Mayer, L.M., Rossi, P.M., 1982. Specific surface area in coastal sediments: relationships with other textural factors. Mar. Geol. 45, 241–252. McKeague, J.A., Day, J.H., 1966. Dithionite and oxalate extractable Fe and Al as aids in differentiating various classes of soils. Can. J. Soil Sci. 46, 13–22. Mikutta, R., Kleber, M., Kaiser, K., Jahn, R., 2005. Review: organic matter removal from soils using hydrogen peroxide, sodium hypochlorite and disodium peroxidisulfate. Soil Sci. Soc. Am. J. 69, 120–135. Mongeon, A., Aherne, J., Watmough, S.A., 2010. Steady-state critical loads of acidity for forest soils in the Georgia Basin, British Columbia. J. Limnol. 69, 193–200. Moore, R.D., Trubilowicz, J.W., Buttle, J.M., 2012. Prediction of streamflow regime and annual runoff for ungauged basins using a distributed monthly water balance model. J. Am. Water Resour. As 48, 32–42. Ouimet, R., Duchesne, L., 2005. Base cation mineral weathering and total release rates from soils in three calibrated forest watersheds on the Canadian Boreal Shield. Can. J. Soil Sci. 85, 245–260. Ouimet, R., Duchesne, L., Houle, D., Arp, P.A., 2001. Critical loads and exceedances of acid deposition and associated forest growth in the northern hardwood and boreal coniferous forests in Quebec, Canada. Water Air Soil Pollut. 1, 119–134. Ouimet, R., Arp, P.A., Watmough, S.A., Aherne, J., DeMerchant, I., 2006. Determination and mapping critical loads of acidity and exceedances for upland forest soils in eastern Canada. Water Air Soil Pollut. 172, 57–66. Petersen, L.W., Moldrup, P., Jacobsen, O.H., Rolston, D.E., 1996. Relations between specific surface area and soil physical and chemical properties. Soil Sci. 161, 9–21. Phelan, J., Belyazid, S., Kurz, D., Guthrie, S., Cajka, J., Sverdrup, H., Waite, R., 2014. Estimation of the base cation weathering rates with the PROFILE model to determine critical loads of acidity for forested ecosystems in Pennsylvania, USA: pilot application of a potential national methodology. Water Air Soil Pollut. 225, 2109.

9

Reid, C.R., Watmough, S.A., 2016. Spatial patterns, trends and the potential long-term impacts of tree harvesting on lake calcium levels in the Muskoka River watershed, Ontario, Canada. Can. J. Fish. Aquat. Sci. 73, 382–393. Santamarina, J.C., Klein, K.A., Wang, Y.H., Prencke, E., 2002. Specific surface: determination and relevance. Can. Geotech. J. 39, 233–241. Soil Classification Working Group, 1998. The Canadian System of Soil Classification. third ed. Agriculture and Agri-Food Canada. Publication 1646. p. 187. Sverdrup, H., Warfvinge, P., 1988. Weathering of primary silicate minerals in the natural soil environment in relation to a chemical weathering model. Water Air Soil Pollut. 38, 387–408. Sverdrup, H., Warfvinge, P., 1992. Critical loads for forest soils in Nordic countries. Ambio. 21, 348–355. Sverdrup, H., Warfvinge, P., 1995. Critical loads of acidity for Swedish forest ecosystems. Ecol. Bull. 44, 75–89. Team, RStudio, 2015. RStudio: Integrated Development for R. RStudio, Inc., Boston, MA Accessible at:. http://www.rstudio.com/. UNECE (United Nations Economic Commission for Europe), 2004. Manual on Methodologies and Criteria for Modelling and Mapping Critical Loads and Levels and Air Pollution Effects, Risks and Trends. UNECE Convention on Long-range Transboundary Air Pollution Accessible at:. www.icpmapping.org/Mapping_Manual. UNECE (United Nations Economic Commission for Europe), 2016. Towards Cleaner Air Scientific Assessment Report 2016: North America. UNECE Convention on Longrange Transboundary Air Pollution Accessible at:. http://www.unece.org/fileadmin/ DAM/env/documents/2016/AIR/Publications/LRTAP_Assessment_Report_-_North_ America.pdf. Warfvinge, P., Sverdrup, H., 1992. Calculating critical loads of acid deposition with PROFILE- a steady-state soil chemistry model. Water Air Soil Pollut. 63, 119–143. White, A.F., Brantley, S.L., 2003. The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and the field? Chem. Geol. 202, 479–506. White, A.F., Blum, A.E., Schulz, M.S., Bullen, T.D., Harden, J.W., Petersen, M.L., 1996. Chemical weathering rates of a soil chronosequence on granitic alluvium: quantification of mineralogical and surface area changes and calculation of primary silicate reaction rates. Geochem. Cosmoch. Acta (14), 2533–2550. Whitfield, C.J., Reid, C., 2013. Predicting surface area of coarse-particle sized soils: implications for weathering rates. Can. J. Soil Sci. 93, 621–630. Whitfield, C.J., Watmough, S.A., Aherne, J., Dillon, P.J., 2006. A comparison of weathering rates for acid sensitive catchments in Nova Scotia, Canada and their impact on critical load calculations. Geoderma. 136, 899–911. Whitfield, C.J., Phelan, J.N., Buckley, J., Clark, C.M., Guthrie, S., Lynch, J.A., 2018. Estimating base cation weathering rates in the USA: challenges of uncertain soil mineralogy and specific surface area with applications of the PROFILE model. Water Air Soil Pollut. 229, 61. Williston, P., Aherne, J., Watmough, S., Marmorek, D., Hall, A., de la Cueva Bueno, P., Murray, C., Henolson, A., Laurence, J.A., 2016. Critical levels and loads and the regulation of industrial emissions in Northwest British Columbia, Canada. Atmos. Environ. 146, 311–323.