Increased carbon capture by a silicate-treated forested watershed affected by acid deposition

Meeting internationally agreed-upon climate targets requires Carbon Dioxide Removal (CDR) strategies coupled with an urgent phase-down of fossil fuel emissions. However, the efficacy and wider impacts of CDR are poorly understood. Enhanced rock weathering (ERW) is a land-based CDR strategy requiring large-scale field trials. Here we show that a low 15 3.44 t ha wollastonite treatment in an 11.8-ha acid-rain-impacted forested watershed in New Hampshire, USA led to cumulative carbon capture by carbonic acid weathering of 0.025–0.13 t CO2 ha over 15 years. Despite a 0.8–2.4 t CO2 ha logistical carbon penalty from mining, grinding, transportation and spreading, by 2015 weathering together with increased forest productivity led to net CDR of 8.5–11.5 t CO2 ha. Our results demonstrate that ERW may be an effective, scalable CDR strategy for acid-impacted forests but at large-scale requires sustainable sources of silicate rock dust. 20


Introduction
The Intergovernmental Panel on Climate Change (IPCC) (Rogelj et al., 2018) Special Report on global warming indicates large-scale deployment of Carbon Dioxide Removal (CDR) technologies will be required to avoid warming in excess of 1.5 °C by the end of this century. Land-based CDR strategies include enhanced rock weathering (ERW), which aims to accelerate the natural geological process of carbon sequestration by amending soils with crushed reactive calcium (Ca) and magnesium 25 (Mg)-bearing rocks such as basalt (The Royal Society and The Royal Academy of Engineering, 2018; Hartmann et al., 2013).
Forests represent potential large-scale deployment opportunities where rock amendments may provide a range of benefits, including amelioration of soil acidification and provisioning of inorganic plant-nutrients to cation-depleted soils (Hartmann et al., 2013;Beerling et al., 2018). Although ERW has not yet been demonstrated as a CDR technique at the catchment scale, a forested watershed experiment at the Hubbard Brook Experimental Forest (HBEF,43° 56'N, in the White 30 Mountains of New Hampshire, USA provides an unusual opportunity for assessing proof-of-concept in this priority research area. The HBEF watershed experiment, designed to restore soil calcium following decades of leaching by acid rain, involved application of a finely ground rapidly-weathered calcium silicate mineral wollastonite (CaSiO3; 3.44 t ha -1 ) on 19 October 1999 to an 11.8-ha forested watershed (SI Appendix) Peters et al., 2004;Shao et al., 2016). Unlike the 35 carbonate minerals (e.g., CaCO3) commonly applied to acidified soils (Lundström et al., 2003), wollastonite does not release CO2 when weathered (Supplementary Information) so is much better suited for CDR (Hartmann et al., 2013). It also has dissolution kinetics comparable to or faster than other calcium-rich silicate minerals such as labradorite found in basalt (Brantley et al., 2008). Thus, the HBEF experiment provides a timely and unparalleled opportunity for investigating the longterm (15 years) effects of ERW on CDR potential via forest and stream water chemistry responses. 40 In the case of ERW with wollastonite, CDR follows as Ca cations (Ca 2+ ) liberated by weathering consume atmospheric CO2 through the formation of bicarbonate (HCO3 -) by charge balance, as described by the following reaction: However, forests in the northeastern USA have experienced acid deposition (Likens and Bailey, 2014), changes in nitrogen cycling (Goodale and Aber, 2001;McLauchlan et al., 2007) and increases in dissolved organic carbon (DOC) fluxes 45 (Cawley et al., 2014) that may affect CO2 removal efficiency by ERW processes. In particular, CO2 consumption as measured by bicarbonate production may be diminished if sulphate (SO4 2-), nitrate (NO3 -), or naturally-occurring organic acid anions (Fakhraei and Driscoll, 2015) (H2A -) in DOC intervene to inhibit the following mineral weathering reactions. For example: CaSiO3 + H2O + H2SO4  Ca 2+ + H4SiO4 + SO4 2-, CaSiO3 + H2O + 2HNO3  Ca 2+ + H4SiO4 + 2NO3 -, CaSiO3 + H2O + 2H3A  Ca 2+ + H4SiO4 + 2H2A -, These environmental effects on stream-water chemistry are well documented at the HBEF (Cawley et al., 2014;Likens and Bailey, 2014;Rosi-Marshall et al., 2016;McLauchlan et al., 2007), and may be exacerbated under future climate change (Sebestyen et al., 2009;Campbell et al., 2009).
Here we exploit the experimental design and long-term monitoring of streamwater chemistry, trees, and soils, for two 55 small forested HBEF watersheds to evaluate the effects of the wollastonite treatment in 1999 on catchment CO2 consumption via inorganic and organic pathways. Further, we examine how biogeochemical perturbations in S, N, and organic carbon cycling affect catchment inorganic CO2 consumption. We consider the forest response, the carbon cost for ERW deployment (mining, grinding, transportation and application), and the net greenhouse gas balance for the treatment. Finally, we provide an initial assessment of the net CDR potential of silicate treatments deployed over larger areas of acidified forest in the 60 northeastern United States.

Site description
The HBEF has a temperate climate with ~1400 mm mean annual precipitation of which up to one third falls as snow (Campbell 65 et al., 2007). The mean temperatures in January and July are -9 °C and 18 °C respectively, and the period from mid-May to mid-September comprises the growing season (Campbell et al., 2007). There are six small southeast-facing watersheds in the HBEF with 20%-30% slopes , including a biogeochemical reference (watershed W6, 13.2 ha, 545-791m asl) and one which received the silicate treatment (watershed W1, 11.8 ha, 488-747m asl). Carbonate and evaporite minerals are in very low abundance (<1% calcite in the crystalline rocks and glacial deposits) in these silicate-mineral 70 dominated watersheds (Johnson et al., 1981). Well-drained Typic Haplorthod soils with pH<4.5 and mean depth 0.6m formed from relatively impermeable glacial till, which restricts water flow and protects the underlying schist bedrock from weathering.
Fagus grandifolia, Betula allegheniensis and Acer saccharum are the dominant trees in this Northern Hardwood forest, while Betula papyrifera, Abies balsamea and Picea rubens are common at the highest elevations where soils tend to be shallow and wetter (Cho et al., 2012). A. saccharum and P. rubens are both calcium-sensitive, but soil calcium-bearing minerals are 80 less available to A. saccharum (Blum et al., 2002) and total bioavailable calcium content decreases with elevation (Cho et al., 2012). This silicate-addition experiment was designed to replace bioavailable calcium which had been stripped from the soils by decades of acid deposition.

Treatment description
On 19 and 21 October 1999, W1 was treated with 344 g/m 2 of pelletized wollastonite (CaSiO3) by a GPS-equipped helicopter 85 with a motorized spreader to ensure even deployment across the catchment, including the 1804 m 2 streambed (Peters et al., 2004). Following treatment, the lignin-sulfonate binder forming the pellets dissolved within several days (Peters et al., 2004), and the ground wollastonite itself dissolved rapidly in the upper Oie soil horizon, increasing Oie base saturation from 40% to 78% and raising soil pH from 3.88 to 4.39 within one year . Although the budget of wollastonite-derived calcium (Wo-Ca) has never been closed due to lack of data from vegetation and from deeper soil layers (Shao et al., 2016), it 90 https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License. is thought that uptake by vegetation and retention by soil exchange sites delayed transport of Wo-Ca to lower soil horizons and streamwater for three years .
Using MATLAB (version R2016a) scripts, we wrote PHREEQC input files and determined the inorganic carbon species for each streamwater sample with PHREEQC. Along with a standard database which decouples ammonium and nitrate 100 (Amm.dat, provided with the PHREEQC software), we included the ionization constants for the organic acid triprotic analogue and the constants for Al complexation described for Hubbard Brook streams (Fakhraei and Driscoll, 2015) in our PHREEQC simulations. These are: pKa1=2. 02,pKa2=6.63,pKa3=7.30,pKAl1=4.07,pKAl2=7.37,pKAl3=6.65, and site density m=0.064 mol sites mol C -1 . Our organic acid concentrations are the product of the corresponding site density of reactions and the measured dissolved organic carbon concentration (Fakhraei and Driscoll, 2015); these were PHREEQC inputs along with total 105 monomeric Al and major ion concentrations from the longitudinal datasets.
Spectator ions (Cland NH4 + ) were adjusted to achieve charge balance given the measured pH for the treated and reference watersheds. Clwas only adjusted when charge balance was not achieved using NH4 + alone. This was deemed to be the case when PHREEQC failed to converge or when the percent error exceeded 5%. We used original rather than adjusted rainwater Cl to calculate the contribution of rainwater to streamwater chemistry (described below). These adjusted ions were 110 then held constant for our modelled scenarios, while pH was allowed to vary.
Exploratory PHREEQC tests (charge-balancing on DIC) either with or without organic acids suggest that the acids depress total DIC, HCO3and also the saturation state of gaseous CO2. Similar variability in the saturation is also observed when DIC values from partially degassed samples from the streams are used as input. We chose minimum and maximum values of 1100 and 1700 ppm, or ~3 and 4.6 × 368, the mean value of Mauna Loa pCO2 for 1985-2012. These values correspond to 115 log10(pCO2(g)) = -2.87±0.09 SD derived from a prior analysis of this variability for the same time range (Fakhraei and Driscoll, 2015).
temperatures. These temperatures were used in our PHREEQC modelling, with equilibrium constants for the DIC species as functions of temperature. Only samples measured closest to the weirs and with a valid pH were processed with PHREEQC.

Catchment CO2 consumption
Annual total watershed CO2 consumption (Eq. 1) was calculated as the product of streamwater flow and [HCO3 -]stream, corrected for the HCO3contribution of rainwater (αrain,HCO3, see below), as given by: where [HCO3 -]stream is given in mol kgw -1 and flow is the "runoff" in mm time -1 . Calculated [HCO3 -]stream and annual CO2 consumption for the treated and reference watersheds (Eqs. 5 and 13) comprise our baseline simulations and represent a primary test of hypothesized increased carbon capture resulting from weathering of the applied silicate.
Bicarbonate-derived CO2 consumption (Eq. 5) is the most conservative approach to estimating net carbon fluxes related 140 to ERW. For natural freshwaters in equilibrium with the atmosphere, this entails a titration for total alkalinity with a possible correction for the concentration of organic acid anions (Köhler et al., 2000). However, another widely used (Jacobson and Blum, 2003) measure of CO2 consumption is derived assuming charge-balance of base cations (Ca 2+ , Mg 2+ , K + and Na + ) by where streamwater cation equivalents are corrected for contributions from rain/snow precipitation (Methods) and sulphuric acid weathering (Chetelat et al., 2008)  the applied minerals may occur, but we consider that the cations released from the applied minerals comprise the most unambiguous treatment effect. The charge associated with wollastonite-derived Ca 2+ (Wo-Ca) determines the CO2 150 consumption associated with the HBEF wollastonite treatment: where X is the fraction of Ca 2+ from weathered wollastonite. Equations 8 and 9 are optimistic measures of CO2 consumption because they ignore both the weathering agent and streamwater inorganic carbon, assuming charge-balance of cations by carbonate and bicarbonate ions in the oceans. Eqs. 7 and 9, together with our flux calculations accounting for sparsity of 155 concentration data compared to daily flow data (Methods), should help avoid major uncertainties in catchment-scale CO2 consumption calculations: the provenance of the cations and variations in concentration and discharge (Moon et al., 2014).

Fraction of calcium derived from wollastonite
We applied a two-component mixing model that , (10)  160 where pre-app and post-app refer to pre-application and post-application streamwater concentrations and Wo refers to wollastonite. The Sr data  have been extended through 2015. See SI Appendix for further discussion of the use of strontium and its isotopes as tracers of Ca 2+ provenance.

Contributions of rain/snow precipitation to streamwater chemistry
We estimated the contribution of rain/snow (Likens, 2016b, a) relative to all other sources, using a previously published mixing 165 model (Négrel et al., 1993). We assume all Clin the water is from rain/snow, noting that this common treatment of Cl as an unreactive tracer is not always justified (Lovett et al., 2005). The contribution of precipitation to the streamwater (αrain) is generally set using Na and Cl, which are less affected by nutrient cycling and adsorption than other major ions (Négrel et al., 1993): To account for attenuation of the rain/snow precipitation leaching through the soil, Cl/Na and HCO3/Na at any given time (t) are means from the previous three months. The contribution of rain/snow to other ions such as HCO3in the streamwater can be estimated as follows: rain,HCO3,t = rain,Na,t × 175 https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License.

Flux calculations
To ensure that fluxes from our two watersheds were comparable and to correct for the sparsity of solute measurements compared to flow measurements, we created rolling annual flow-adjusted fluxes using Method 5 of Littlewood et al (Littlewood et al., 1998) at five evenly-spaced points each year: 180 where Qi is the measured instantaneous stream flow, Ci is the concentration for sample i, M is the number of streamwater chemistry samples in the year (usually 12), Qk is the k th flow measurement, and N is the number of flow measurements. In our case, daily flow measurements (Campbell, 2015) and ~monthly streamwater samples (Driscoll, 2016a, b) were available. 185 Therefore, the mean concentration for the preceding twelve months is multiplied by the mean flow for the same period, suitably scaled to get the total annual flux. Without sub-daily timestamps for the longitudinal streamwater chemistry data, we used daily total flows rather than instantaneous flows. Tests suggested that there was little difference between using mean daily instantaneous flows and the mean daily total flows.

Carbon sequestration in wood
Battles et al  provided mean wood production over two five-year periods for the treated and reference watersheds. We considered the difference (treated-reference) to be an estimate of the treatment effect on potentially long-term biomass carbon sequestration. Assuming 46.5% of the woody biomass is carbon (Martin et al., 2018), our calculated cumulative additional C sequestration in the treated watershed over ten years [given by (5×0.78 + 5×0.29) × 0.465 × 100 / 12] 195 was 20.7 mol C m -2 (9.1 t CO2 ha -1 ).

Greenhouse gas emissions from soils
Measurements (Groffman, 2016) were taken at four elevations in the treated watershed and at points just west of the reference watershed starting in 2002. Gas samples were collected from chambers placed on three permanent PVC rings at each of these eight sites (Groffman, 2016). The data were not normally distributed so were analyzed with Kruskal-Wallis tests at the 0.05 200 significance level; however, tests with one-way ANOVA produced the same overall results. All analyses were done in Matlab R2016a.
Cumulative curves for each of the 24 chambers were generated by matching the dates of the measurements, excluding points which were missing data for any chamber and allowing up to a week's discrepancy between catchments. Nearly all discrepancies were within one day. Assuming diurnal variation was minor compared to seasonal variation, each datum (g C 205 https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License. m -2 hour -1 ) was multiplied by 24 hours and by 30 days to get gC m -2 month -1 . There was no extrapolation to fill gaps in the dataset; curves shown in Fig. 3 are internally consistent but not comparable to other datasets. We were particularly interested in the elevation-specific responses, as the different elevations have distinct tree species compositions and below-ground responses to the wollastonite treatment .
The HBEF experimental watersheds are divided into 25×25m plots on slope-corrected grids. Vegetation has been 210 surveyed four times since the late 1990s and assigned a zone designation in each plot Driscoll Jr et al., 2015;Battles et al., 2015b, a) (Supplementary Fig. S9). To estimate the respiration savings over the whole watershed, we added the areas of individual plots which were assigned to our four vegetation types (Low, Mid and High hardwoods, and Spruce-Fir). Because there were seven vegetation types in the datasets, we compared all types with pairwise Kruskal-Wallis tests at the 0.05 significance level using the basal area data for the six dominant tree species. Kruskal-Wallis tests were 215 appropriate because the data, and therefore the differences from the means (residuals), were not normally distributed. These tests suggested that the "extra" vegetation types ("Birch/Fern Glade", and "Poor Hardwoods" at High and Mid elevations) could be combined with Spruce-Fir, High and Mid Hardwoods respectively. Watershed fractions for our combined forest types were 0.155 for SpruceFir, 0.16 for High Hardwoods, 0.415 for Mid Hardwoods, and 0.27 for Low Hardwoods. When creating our composite treatment effects for the entire watershed, we considered a treatment effect to be present only where our 220 statistical analyses suggested significantly different fluxes.

Logistical carbon emissions costs
We used the 1999 upstate New York CO2 emission factor for electricity generation from oil (United States Environmental Protection Agency, 1999) (0.9 Mg CO2 MWh -1 ), and rearranged Equation 28 of Stamboliadis (Stamboliadis et al., 2009) where the specific surface area s (1600 m 2 kg -1 for our treatment) is related to the specific potential energy ep of the material (kJ kg -1 ), with theoretical parameters (Stamboliadis et al., 2009) α=139 m 2 kJ -1 and μ=0.469 (dimensionless). We convert this potential energy to MWh t -1 Qz (3600 seconds per hour and 1000 kWh MWh -1 ). The equation was derived for quartz (Qz) which has hardness 7. Because wollastonite hardness is in the range 5-5.5, this equation may overestimate the energy needed to grind the wollastonite. 230 The main energy source in Allerton will have been coal, and the 1999 Illinois emissions factor (United States Environmental Protection Agency, 1999) is 1.1 Mg CO2 MWh -1 . The monetary cost is USD0.041 kWh -1 for pelletization of limestone fines and USD0.85 t -1 product, so we estimate 20.73 kWh t -1 product.
HBEF experiment are provided in Table 2 with calculation details given in Table 3. The Matlab script used for these calculations is available on request. Note: t refers to megagrams, not US short tons.

Greenhouse gas budget for a treatment
The success of any treatment for climate change mitigation is determined by the net greenhouse gas (CO2 equivalent) fluxes 240 prior to and following treatment, at the treatment site and downstream. The GHG balance associated with any treatment, including ERW, is given by where the individual fluxes are described and values given in Table 2. Here, ΔTOC and ΔNO3N2O are penalties because these lead to CO2 and N2O emissions downstream. 245 Changes in net ecosystem productivity due to the treatment ΔNEP can be represented by where SOC is soil organic carbon, SRESP is soil (root+heterotrophic) respiration and ARESP is aboveground (wood+canopy) respiration. For the HBEF wollastonite experiment, we lack belowground biomass and aboveground respiration. Wood is a longer-term carbon sink than leaves or twigs so we have chosen to let this represent our biomass increment. Eq. (15) neglects 250 ecosystem disturbances including fire, and possible carbonate mineral precipitation in soils. There is no evidence for the latter at the HBEF.
With the exception of CO2 consumption, treatment effects were calculated as the difference (treated-reference) between the two watersheds. We used a range of emissions factors for N2O to estimate the penalty associated with nitrate export (ΔNO3N2O); low: 0.0017 kgN2O-N kg -1 DIN (Hu et al., 2016) and high: 0.0075 kgN2O-N kg -1 DIN (De Klein et al., 2006), 255 where DIN is dissolved inorganic nitrogen dominated by nitrate. This N2O was then converted to CO2e (CO2 equivalents in terms of cumulative radiative forcing) given the 100-year time horizon global warming potential (Pachauri et al., 2014) (GWP100) for N2O: 265 gCO2e g -1 N2O. Likewise, ΔCH4 was converted to CO2e (CO2 equivalents in terms of cumulative radiative forcing) given GWP100 for CH4: 28 gCO2e g -1 CH4.

Wollastonite treatment increased streamwater CO2 export
We first consider the time-series of streamwater changes in Ca 2+ concentrations in the treated ([Ca]T,stream) and reference ([Ca]R,stream) watersheds. Immediately after treatment, [Ca]T,stream increased from <30 μmol L -1 to ~60 μmol L -1 , and then slowly declined over the next decade, remaining persistently above [Ca]R,stream for 15 years (Fig. 1a). The initial post-treatment peak represents dissolution of wollastonite within the stream (Peters et al., 2004) and release of calcium from hyporheic exchange 265 during the first few years (Shao et al., 2016;Nezat et al., 2010). Retention of Ca 2+ ions liberated by wollastonite dissolution We derived the annual export of Ca 2+ from the treated and reference watersheds as the product of mean annual flowadjusted Ca 2+ streamwater concentrations and annual flow (Fig. 1b) (Methods). After accounting for variations in flow, increased streamwater Ca 2+ concentrations in the treated watershed are translated into a 2-fold increase in total Ca 2+ export relative to the reference watershed that was maintained for 15 years until 2015 through this analysis period. Overall, the 275 wollastonite treatment resulted in a sharp spike in calculated CO2 consumption (Wo-CO2,Ca) that decreased but remained elevated as a result of the treatment (Fig. 1c).
Temporal patterns in modelled streamwater bicarbonate concentration in both treated and reference watersheds (Fig.   1d), and the corresponding total annual CO2 consumption (CO2,HCO3) (Fig. 1e) and CO2 consumption resulting from treatment (Wo-CO2,HCO3) (Fig. 1f), largely mirror changes in streamwater Ca 2+ concentrations but are modified by the supply and loss 280 of anions. Calculated flow-adjusted CO2 consumption (Fig. 1e) peaked 2-3 years post-treatment with a broader peak in CO2 consumption evident in 2007-2012 corresponding to declining legacy effects of acid rain until transient NO3peaks appeared 2012-2015. The Wo-CO2,HCO3 shows a pattern that mirrors the Wo-CO2,Ca but is generally 5 times lower (Fig. 1c,f).

Sulphuric, nitric and organic acids reduce CDR
We next undertook sensitivity analyses to investigate the effects of acid deposition, increased NO3and organic acid export 285 from the treated watershed on bicarbonate concentrations and resulting CO2 consumption (Fig. 2). In a 'Low SO4' scenario ( Fig. 2a-c), we sought to understand the effects of acid deposition by replacing the mean monthly time-series of streamwater and rainwater SO4 2for the treated watershed with a new time-series (purple curve, Fig. 2a) created by repeating the post-2010 datasets, which reflect diminished acid deposition following emission controls from the US Clean Air Act (Likens and Bailey, 2014). Removing acid rain effects in this manner dramatically increased the calculated bicarbonate concentrations and total 290 annual CO2 consumption (CO2,HCO3), increasing the initial spikes resulting from the wollastonite treatment in both by at least four-fold (purple curves, Fig. 2 b,c). An additional legacy of acidification in North American forests (Harrison et al., 1989) is SO4 2retention on soil clay mineral Fe and Al oxides (Fuller et al., 1987), which were subsequently released by increased soil pH following wollastonite weathering (Shao et al., 2016;Fakhraei et al., 2016). To assess the effect of this legacy SO4 2-, we ran simulations substituting the lower streamwater SO4 2concentrations from R in T (T REF, green curves, Fig. 2b,c). Results 295 suggest that legacy SO4 2accounts for over half of the total acid deposition effect on increased [HCO3 -]stream and CO2 consumption in the simulations.
In the 'Ref NO3' scenario ( Fig. 2 d-f), seasonal spikes in streamwater export of NO3recorded from the treated watershed between 2012 and 2015 were removed by substituting the reference watershed streamwater NO3concentration measurements https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License. lacking these spikes. This manipulation markedly increased modelled bicarbonate (Fig. 2e) and mean annual CO2 consumption 300 (Fig. 2f). To quantify the effects of organic acids on bicarbonate production in the treated watershed, we ran "+OA" and "-OA" simulations, i.e., with and without accounting for organic acids, respectively (Fig. 2 g-i). Results showed that removing OA from our simulations also increased modelled streamwater bicarbonate concentration (Fig. 2h), and resulting CO2 consumption (Fig. 2i), in the treated watershed.

Effects of increasing wollastonite treatment 305
Because the HBEF application rate (3.44 t ha -1 ) is smaller than the 10-50 t ha -1 suggested for ERW strategies (Strefler et al., 2018;Beerling et al., 2018), we simulated the possible effects of a ten-fold increase in the streamwater Ca 2+ concentrations on bicarbonate production (Fig. 3a) and CO2 consumption (Fig. 3b). In this initial assessment, we assume streamwater responses are directly proportional to wollastonite application rate, i.e., 34.4 t ha -1 , and that all other variables remained unchanged.
Results show that after 15 years, cumulative Wo-CO2,HCO3 is 73% of Wo-CO2,Ca (Fig. 3c), as opposed to less than 20% for the 310 actual rate of 3.44 t ha -1 ( Table 2). These results suggest that at higher application rates of wollastonite, the details of the CO2 consumption calculations become less important.

Amplification of organic carbon sequestration by wollastonite treatment
In reversing long-term Ca 2+ depletion of soils, the silicate rock treatment significantly increased forest growth and wood production between 2-12 years post-treatment relative to the reference watershed . This forest response 315 increased total carbon sequestration by 20.7 mol C m -2 or 9.1 t CO2 ha -1 during those ten years as a result of the treatment (Methods).
Changes in greenhouse gas (GHG) emissions from soils represent a further route to affecting the climate mitigation potential of the wollastonite treatment. Despite a rapid increase of one pH unit in the upper organic soil horizon (Oie), soil respiration CO2 fluxes showed no significant difference between watersheds during the first three years after treatment 320 . However, our analysis of newly available longer-term datasets indicates that the treatment significantly reduced soil respiration in the high elevation hardwood zone (~660-845m a.s.l.) (χ 2 (1,270)=17.2, P < 0.001), possibly due to reduced fine-root biomass  rather than changes in microbial activity . No significant effects on soil respiration were detected in any of the other HBEF vegetation zones (Fig. 4). The wollastonite treatment increased the soil sink strength for CH4 (χ 2 (1,266)=30.8, P < 0.001) in the low-elevation hardwood zone (482-565m 325 a.s.l.), while it decreased in the high elevation zone (χ 2 (1,268)=22.3, P < 0.001) (SI Appendix, Fig. S8). There were no significant treatment effects on soil N2O fluxes in any vegetation zone (SI Appendix).
Results show that after 15 years, cumulative Wo-CO2,HCO3 is 73% of Wo-CO2,Ca (Fig. 3c), as opposed to less than 20% for the actual rate of 3.44 t ha -1 ( Table 2). These results suggest that at higher application rates of wollastonite, the details of the CO2 consumption calculations become less important. 335

Potential for deployment at larger scales
The HBEF forests are representative of a major area of eastern North America receiving acid deposition since the 1950s (Likens and Bailey, 2014) which may be suitable for remediation and carbon capture via ERW treatment with a silicate rock or mineral. For example, the Appalachian and Laurentian-Acadian Northern Hardwood Forests (NHWF) covering a combined area of 0.137 Mkm 2 in the United States (Ferree and Anderson, 2013) have the same dominant hardwood trees as the HBEF 340 experimental watersheds (Fagus grandifolia, Betula allegheniensis and Acer saccharum). Acid deposition exceeded "critical loads" likely to harm ecosystems in almost 9000 ha of New Hampshire's Acer saccharum stands (NHAs) (Schaberg et al., 2010). These forests might be expected to respond similarly to a wollastonite treatment. The acid-sensitive trees Acer saccharum and Picea rubens are also widely distributed along the high elevation acid sensitive regions of the Appalachian Mountains which have already been impacted by acid deposition (Lawrence et al., 2015). We define this as a 40-km corridor 345 along the Appalachian Mountains comprising 0.14 Mkm 2 and overlapping with the High Allegheny Plateau Ecoregion (HAL) where Acer saccharum is declining above ~550 m a.s.l. ) (0.07 Mkm 2 ).
We examined the potential CO2 consumption for a range of wollastonite application rates encompassing those suggested for ERW strategies (Strefler et al., 2018;Beerling et al., 2018) (Fig. 6). In this analysis, we adjusted mean (2003-2012) Wo-CO2,ca for the actual 3.4 4 t ha -1 treatment (~0.2 mol C m -2 yr -1 ) proportionally for 10-50 t ha -1 treatments. We assume logistical 350 carbon penalties are minimised and balanced by forest biomass carbon sequestration responses to treatment. This analysis suggests net CDR potential of 0.3-1.7 Mt CO2 yr -1 along the Appalachian corridor, which is 2-12% of New Hampshire state emissions (13.8 Mt CO2) in 2016 (Energy Information Administration, 2019). However, world wollastonite reserves (Curry, 2019) (≥0.1 Pg) are insufficient to treat large areas of eastern North America at rates of 10-50 t ha -1 , highlighting the requirement for alternative sustainable sources of silicate materials. 355

Discussion
Our analyses of wollastonite application at the HBEF provide a unique long-term (15 year) perspective on the whole watershed carbon cycle responses and net CDR by accounting for the associated CO2 costs of logistical operations. By 2015, net CDR amounted to 8.5-11.5 t CO2 ha -1 at a low rate of wollastonite application, with increased carbon sequestration into forest biomass playing the dominant role. We estimate that if the HBEF application rates were increased ten-fold, net CDR would 360 increase by 8%, assuming 400-km transport distances and no change in forest responses. Amplification of organic carbon https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License. capture may therefore represent a major CDR benefit of ERW when applied to forested lands affected by acid rain. Forest management practices, disturbance regimes and the ultimate fate of any harvested wood are also important in determining the storage lifetime of the sequestered carbon. Our results highlight the need to carefully monitor the net carbon balance of forested ecosystems in response to a silicate treatment, including wood and canopy respiration (Fahey et al., 2005) (Methods). This 365 challenging goal might best be achieved with fully instrumented eddy covariance plots, although the HBEF topography is not well suited for this approach (Fahey et al., 2005).
Inorganic CO2 consumption calculated based on streamwater bicarbonate fluxes approximately doubled in the treated watershed relative to the reference watershed 15 years post-treatment (0.028 and 0.016 tCO2 ha -1 , respectively) ( Table 1). The presence of SO4 2-, NO3and organic acid anions lowered the efficiency of CO2 consumption by alkalinity generation, with acid 370 deposition having the single largest calculated effect ( Table 1). The cause of increased NO3export from the treated watershed is not as yet understood (Rosi-Marshall et al., 2016). If it proves a general feature of terrestrial ecosystem responses to silicate mineral treatment, this could affect the efficiency of carbon capture via bicarbonate export. Overall, we suggest that continued recovery of eastern North American and European forests and soils from acid deposition creates conditions beneficial to watershed health, carbonic acid-driven weathering and inorganic carbon export following application of crushed silicate 375 minerals.
In Asia, acid rain is an ongoing problem with an estimated 28% of Chinese land area (~2.7 Mkm 2 ) receiving potentially damaging S deposition in 2005 (Zhao et al., 2009), and critical loads were exceeded in ~0.36 Mkm 2 of the European Economic Area (EEA) in 1999 , approximately double the affected area of US Northern Hardwood Forests (Fig. 6). Fig. 6 suggests that a single 30t Wo ha -1 treatment over 0.14 Mkm 2 (Appalachian Trail corridor) could, in principle, sequester 380 ~1 MtCO2 y -1 or 15 MtCO2 over 15 years via wollastonite-derived Ca export in streamwater alone. Adding the Chinese and European acidified areas could potentially sequester 0.34 GtCO2, approximately 0.2-0.7% of the ~50-150 Gt CDR required by 2050 to avoid warming in excess of 1.5° (Rogelj et al., 2018). Inclusion of biomass and soil responses increases CDR contributions from ERW on acidified forests, but these will still be modest. Assuming no further forest responses beyond the 15-year HBEF timeframe, we report a GHG balance of ~10 tCO2e ha -1 . This translates to 1 GtCO2e Mkm -2 suggesting 3.2 385 GtCO2e over 15 years for the Appalachian Trail, the EEA and China combined, or 2-6% of global required CDR as described above.
It is uncertain whether other acidified forest ecosystems would respond similarly to the HBEF Acer saccharum forests in New Hampshire. Many Chinese soils (Duan et al., 2016), as well as old deep soils in areas such as the Virginian Blue Ridge Mountains and the German Harz and Fetchel Mountains (Garmo et al., 2014) have high SO4 2sorption capacity. These soils 390 may retain substantially more SO4 2than the HBEF soils, with potential for prolonged SO4 2flushing following ERW treatment and lower bicarbonate production. Liming studies suggest a range of other effects, some of which may also occur with silicate treatments. Liming increases nitrate export, migration of heavy metals and acidity to deeper soil, and fine root production in topsoils leading to frost damage (Huettl and Zoettl, 1993). https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License.
Many forests have been limed with carbonate minerals such as calcite and dolomite to mitigate acidification in the past. 395 Dolomite has also helped reverse Mg deficiency in conifers (Huettl and Zoettl, 1993). Liming generally improves water quality, although it also forms mixing zones with high-molecular-weight Al complexes toxic to fish (Teien et al., 2006). With silicate treatments, nontoxic hydroxyaluminosilicates form instead (Teien et al., 2006). Unfortunately, carbonates are contraindicated for CDR on acid soils because they can be a net source of CO2 in the presence of strong acids (Hamilton et al., 2007). Treatments of European and North American acidified forests with calcite (1-18 t ha -1 CaCO3) or dolomite (2-8.7 t 400 ha -1 CaMg(CO3)2) have, in general, resulted in increased DOC export and soil respiration without increasing tree growth, regardless of forest composition (Lundström et al., 2003). As calcite and dolomite are 44% and 48% CO2 by weight, these treatments will have released 0.44-7.9 and 0.96-4.54 t CO2 ha -1 respectively when fully dissolved, although dissolution may be slow. Over six years following a 2.9 t dolomite ha -1 treatment (90% 0.2-2.0 mm grains) in a Norwegian coniferous watershed equating to 1.36 t CO2 ha -1 , less than 1% of the dolomite dissolved (Hindar et al., 2003). We estimate that CO2 405 consumption corrected for CO2 release and as measured with dolomite-derived Ca and Mg in streamwater (Dol-CO2,Ca+Mg) averaged 0.02 mol CO2 m -2 yr -1 . CO2 release from carbonate minerals equals Ca and Mg release on a molar basis, so 0.02 mol Dol-CO2 m -2 yr -1 was also either exported in streamwater or lost to the atmosphere. This experiment may have a negative greenhouse-gas balance depending on logistical penalties and soil respiration, as there was no significant treatment effect on tree growth or vitality (Hindar et al., 2003). Ca-sensitive Acer saccharum is present at Woods Lake in New York State, yet 410 tree biomass decreased with no significant differences relative to reference catchments during the 20 years following a 6.89 t Mg-calcite ha -1 application (Melvin et al., 2013), equivalent to 3.07 t CO2 ha -1 given 8% Mg content of the pellets. In contrast to our study and other liming studies, root biomass and soil carbon stocks increased in response to this treatment, although soil respiration was reduced (Melvin et al., 2013). Acer saccharum basal area and crown vigour increased over 23 years in response to 22.4 t dolomitic limestone ha -1 (equivalent to 10.0 t CO2 ha -1 ) on the Allegheny Plateau, although basal area and survival of 415 another dominant canopy species, Prunus serotina, was reduced (Long et al., 2011). Clearly, forest responses to mineral treatments are species-and site-specific.
Although the HBEF experiment used wollastonite, this is not a target mineral for ERW, both because of its limited reserves (Curry, 2019) and high monetary costs (Schlesinger and Amundson, 2018). Recent all-inclusive guide prices of ~700 USD Mg -1 for helicopter deployment of pelletized lime along the Appalachian Mountain corridor are comparable to the price 420 of 694 USD Mg -1 for unpelletized 10-µm wollastonite in 2000 (Virta, 2000). Less expensive materials such as locally-sourced waste fines from mines should be considered if their heavy metal content is low, but the choice of treatment material should be considered together with the vegetation and the native minerals. Application of magnesium-rich materials (e.g. olivine), for example, may help reverse Mg deficiency in Pinus sylvatica and Picea abies as dolomite has done (Huettl and Zoettl, 1993), but some other tree species, such as Acer saccharum, have a higher demand for calcium than for magnesium (Long et 425 al., 2009). The treatment of ecologically sensitive catchments always requires caution as some species, such as Sphagnum mosses and lichens, may respond poorly to treatment (Traaen et al., 1997).

Conclusions
This study identifies key challenges in accounting for CO2 removal and feasibility for larger-scale rollout of ERW. Our two methods for calculating the treatment effect on CO2 consumption depend on strontium isotope data and a site-specific mixing 430 model. At the HBEF, forest responses helped repay the initial carbon treatment penalties, but this depended on the presence of a Ca-sensitive species growing on soils stripped of calcium by decades of acid deposition. Given good access roads, helicopter deployment and pelletization may be unnecessary for application on land undergoing afforestation. Co-deployment of ERW with the CDR techniques favoured by the IPCC (Pachauri et al., 2014), afforestation and bioenergy with carbon dioxide capture and storage (BECCS), may contribute to the twin aims of mitigating the effects of acid deposition on forests and CDR, but 435 site-specific research is required to assess the efficacy and suitability of such strategies.

Code availability
The aqueous geochemistry software PHREEQC software, along with documentation, is freely available from the USGS website (https://www.usgs.gov/software/phreeqc-version-3). MATLAB ® may be purchased from the MathWorks website (https://uk.mathworks.com/products/matlab.html). the scripts used to process the HBEF data are available from the 440 corresponding author, without guarantees that these will run with MATLAB versions other than R2016a.

Data availability
Our data are available from the Long Term Ecological Research (LTER) Network Data Portal. This public repository can be accessed via the Hubbard Brook Ecosystem Study website: https://hubbardbrook.org/d/hubbard-brook-data-catalog See Supplement for a full list of filenames, package IDs, DOIs and access dates.

Competing interests
The authors declare that they have no conflict of interest.  https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License. (a) Calcium and bicarbonate concentrations, along with observed reference calcium. Calcium is charge-balanced by up to two moles of 675 bicarbonate. (b) CO2 consumption due to this higher treatment (Wo-CO2,Ca and Wo-CO2,HCO3) and (c) Cumulative CO2 consumption for both this higher treatment and the actual treatment (Wo-CO2,Ca). We have assumed that a 10-fold higher treatment produces 10-fold higher calcium concentrations with no change in sulphate, nitrate or DOC. Simulated HCO3concentration and Wo-CO2,HCO3 account for the presence of organic acids (+OA) given observed DOC. All CO2 consumption curves (b,c) were flow-normalised and corrected for sparsity of samples (Methods). 680 https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License.

Figure 4: Long-term soil respiration responses to wollastonite treatment at Hubbard Brook Experimental Forest.
Cumulative soil CO2 respiration responses of treated and untreated (a) high elevation hardwoods, (b) high elevation conifers, (c) low elevation hardwoods or (d) mid-elevation hardwoods. Plots show cumulative means ± 1 SE for three chamber measurements at each site and time. Reference data were collected from untreated forests immediately adjacent to the western edge of our reference catchment. P-values from Kruskal-Wallis 685 tests comparing treated and reference raw data (SI Appendix) are shown.
https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License.  (Table 2). These include our calculated treatment effects at the end of 2014 along with the wood production produced over ten years published by Battles et al . The CO2 consumption range is given by Wo-CO2,HCO3 calculated by Eq. (7) and Wo-CO2,Ca calculated by Eq. (9). Nitrate export in streamwater leading to N2O greenhouse gas emissions downstream and a small increase in the soil CH4 sink have been converted to CO2-equivalents (Methods). We found no significant effect for soil N2O emissions.

695
Exported DOC is assumed to be respired downstream. Note this is not a full greenhouse gas balance for the treatment as it neglects potential changes in large fluxes such as aboveground autotrophic respiration. (c) Smaller treatment effects from (b). (d) Elements of the greenhouse gas balance associated with a hypothetical ten-fold higher treatment. Here, we assume no change in forest or streamwater treatment responses other than CO2 consumption. 700 https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License. Figure 6: Projected CO2 consumption following higher-dosage treatments. We considered the possibility of higher-dosage silicate treatments on other northeastern United States higher-altitude forests affected by acid rain, such as Acer saccharum forests in New Hampshire (NHAs), the High Allegheny Plateau Ecoregion (HAL), the Appalachian trail corridor (AT), or Northern Hardwood forests (NHWF) dominated by the same tree species as at Hubbard Brook. Because the world's wollastonite reserves (yellow disks) are insufficient to treat 705 these areas, other calcium-rich silicate minerals would be required. CO2 consumption due to higher dosage (t ha -1 ) is estimated as: (mean observed CO2,Ca between 2004 and 2012) × area × dosage / 3.44 t ha -1 .  3154 km = 1397 km (Gouverneur to Allerton) + 1757 km (Allerton to Woodstock) Transport distance for "local pelletization" calculation (km) 408 km (Gouverneur to Woodstock) https://doi.org/10.5194/bg-2020-288 Preprint. Discussion started: 31 July 2020 c Author(s) 2020. CC BY 4.0 License.