Biogeosciences Rates of consumption of atmospheric CO 2 through the weathering of loess during the next 100 yr of climate change

Quantifying how C fluxes will change in the future is a complex task for models because of the coupling between climate, hydrology, and biogeochemical reactions. Here we investigate how pedogenesis of the Peoria loess, which has been weathering for the last 13 kyr, will respond over the next 100 yr of climate change. Using a cascade of numerical models for climate (ARPEGE), vegetation (CARAIB) and weathering (WITCH), we explore the effect of an increase in CO 2 of 315 ppmv (1950) to 700 ppmv (2100 projection). The increasing CO2 results in an increase in temperature along the entire transect. In contrast, drainage increases slightly for a focus pedon in the south but decreases strongly in the north. These two variables largely determine the behavior of weathering. In addition, although CO 2 production rate increases in the soils in response to global warming, the rate of diffusion back to the atmosphere also increases, maintaining a roughly constant or even decreasing CO 2 concentration in the soil gas phase. Our simulations predict that temperature increasing in the next 100 yr causes the weathering rates of the silicates to increase into the future. In contrast, the weathering rate of dolomite – which consumes most of the CO 2 – decreases in both end members (south and north) of the transect due to its retrograde solubility. We thus infer slower rates of advance of the dolomite reaction front into the subsurface, and faster rates of advance of the silicate reaction front. However, additional simulations for 9 pedons located along the north–south transect show that the dolomite weathering advance rate will increase in the central part of the Mississippi Valley, owing to a maximum in the response of vertical drainage to the ongoing climate change. The carbonate reaction front can be likened to a terrestrial lysocline because it represents a depth interval over which carbonate dissolution rates increase drastically. However, in contrast to the lower pH and shallower lysocline expected in the oceans with increasing atmospheric CO 2, we predict a deeper lysocline in future soils. Furthermore, in the central Mississippi Valley, soil lysocline deepening accelerates but in the south and north the deepening rate slows. This result illustrates the complex behavior of carbonate weathering facing short term global climate change. Predicting the global response of terrestrial weathering to increased atmospheric CO2 and temperature in the future will mostly depend upon our ability to make precise assessments of which areas of the globe increase or decrease in precipitation and soil drainage.


Introduction
Today, continental weathering consumes about 0.3 Gt of atmospheric CO 2 each year (Gaillardet et al., 1999).This value is close to the preindustrial net carbon exchange between the atmosphere and the terrestrial biosphere (0.4 Gt C yr −1 ) (Solomon et al., 2007).Continental weathering is sensitive to land use change (Raymond et al., 2008) and ongoing climate change (Gislason et al., 2009;Beaulieu et al., 2012).Despite its magnitude, apparent sensitivity, and acknowledged importance over geologic timescales, continental weathering is generally neglected in studies dealing with the evolution of the global carbon cycle and climate over century timescales.
Continental weathering is a function of the mineralogy of the exposed surfaces (Amiotte-Suchet et al., 2003) many other parameters such as physical erosion (West et al., 2005), climate (White and Blum, 1995), and vegetation cover (Le Hir et al., 2011;Taylor et al., 2012).This multi-parameter dependence makes it difficult to assess the response of continental weathering to climate change on a global basis.Parametric laws can be used to predict weathering from mean annual temperature and runoff for several average lithological types (Hartmann et al., 2009).These laws were established from field observations.However, as pointed out in Goddéris et al. (2009), such a parametric method only gives a snapshot of the weathering system, and does not capture the decadal to centennial dynamics.In this contribution, we use mechanistic models to predict the response of weathering to increasing CO 2 in the future by focusing on the natural experiment of loess that was deposited approximately 13 kyr ago along the Mississippi valley transect in the USA.We have previously shown that mechanistic weathering models could successfully reproduce mineralogical patterns observed in the loess along the climosequence (Williams et al., 2010;Goddéris et al., 2010).This latter modeling effort lends credibility to our attempt here to predict weathering forward in time.
In addition, loess and loess-derived sediments cover large areas of the continents, especially in the Northern Hemisphere.Ten percent of the continents are presently covered by loessic formations derived from glacial eolian deposits (Rousseau, 2001;Haase et al., 2007).Given their abundance at the Earth's surface, they should contribute significantly to the global weathering flux in the ocean for at least two reasons: (1) loess minerals often have only slightly weathered surfaces, and (2) the specific surface area is high owing to the fine-grained nature of loess.Furthermore, since loess contains highly reactive carbonate phases, loess weathering could potentially respond rapidly to environmental change.Indeed, researchers have suggested that the high surface area of loess may cause fast weathering that can be an important cause for the long-term drawdown of CO 2 (Anderson et al., 1997).Understanding how these loessic formations are weathered might thus be key toward quantifying the future evolution of the global carbon cycle.
In North America, loess is one of the most extensive surficial deposits, especially in the central Great Plains region where it is known as the Peoria loess (Rousseau, 2001;Muhs et al., 2008).Minerals include primary silicates (potassium feldspar, plagioclase feldspar, quartz), secondary clay minerals (dominated by smectites), and carbonate phases (calcite and dolomite) (Muhs et al., 2008).Given this mixed mineralogy, a quantitative description of the evolution of loess weathering as a function of climate change is not an easy task.
Using the WITCH numerical model forced by the GENE-SIS climate model, we previously calculated the weathering rates of the Peoria loess over the last 10 kyr for two sites, one located at the northern end of the Mississippi Valley, and the other located at the southern end (Goddéris et al., 2010).The model starts with the initial mineralogical composition of the aeolian deposits 10 kyr ago.It reproduced the time evolution of the mineralogical composition of the two pedons, simulating mineralogical concentrations close to measured values (Goddéris et al., 2010).Features that were simulated included the deepening of the reaction fronts of dolomite, albite, and potassium feldspar.These results suggest that the calculated weathering rates of the pedons were correctly simulated as a function of climate by the WITCH model.They also show that the present mineralogical composition of the loessic formation is the result of the weathering of the loess under variable climatic conditions.In other words, climate has left its imprint on the loessic profiles.The question we ask in this contribution is: what will be the effect of the climate during the next century on the weathering processes at play in the Peoria loess?We thus start the numerical modeling from the present-day state of the weathering system and project it forward to 2100 for the same two pedons to explore how weathering of Mississippi Valley loess will respond to future climate changes.

Model settings
To project weathering into the future, we ran an atmospheric global circulation model (ARPEGE, Salas y Mélia et al., 2005;Gibelin and Déqué, 2003) under the A1B IPCC carbon emission scenario (Solomon et al., 2007).The output of this simulation (mean air temperature and diurnal temperature range, rainfall, relative air humidity, cloud cover and surface horizontal wind speed) was then used to force a global dynamic vegetation model (CARAIB, Dury et al., 2011).CARAIB predicts the vegetation cover, the biospheric productivity, the below-ground hydrology and CO 2 production along the N-S transect of the Mississippi Valley.These physical variables were then used to force WITCH, a processbased model of weathering reactions (Goddéris et al., 2006;Goddéris et al., 2010; see Appendix A for more details).
The three numerical models have been previously validated.The climate model ARPEGE has been part of the international stretched-grid model intercomparison project (SG-MIP).ARPEGE was able to reproduce the North American climate evolution over 12 yr (1987-1998 period) (Fox-Rabinovitz et al., 2006).The CARAIB biospheric model is a dynamic vegetation model initially developed for the global scale (Warnant et al., 1994).As such it has been used and validated over all continents (Nemry et al., 1996(Nemry et al., , 1999)).Model vegetation distribution maps have also been validated globally (Otto et al., 2002).Here, a globally validated classification of plant functional types (PFT) is used (updated from Franc ¸ois et al., 2011).Finally, the WITCH weathering model has been used to simulate water-mineral interactions in various environments, from polar to tropical settings (Goddéris et al., 2006;Roelandt et al., 2010;Violette et al., 2010;Beaulieu et al., 2011).WITCH was able to reproduce the concentration in base cations and silica in soil solutions and at the outlet of the simulated catchments.
The available ARPEGE climate simulation covers the 1950-2100 period.Instead of running the model cascade only into the future (from 2012 to 2100), we decided to perform the simulation over the whole 1950-2100 period for two reasons: (1) since the weathering response is a priori assumed to be low (specifically for silicate minerals), a longer run will maximize the weathering changes, and (2) the simulation must be started from the less anthropogenically perturbed state, for the purpose of relaxation of the initial condition.

Mineralogical composition of the pedons
In our previous study (Goddéris et al., 2010), we relied on published chemical analyses and X-ray diffraction data for the Peoria loess sampled from pedons along the Mississippi transect (Muhs et al., 2001).Williams et al. (2010) also reported new quantitative X-ray diffraction (XRD) analyses of Muhs' samples from pedons 1 (Clayton County, IL), 13 (Lauderdale County, TN), and 22 (St.Francisville, LA).Those samples were interpreted to contain montmorillonite, illite, plagioclase, potassium feldspar, quartz, and minor kaolinite.In addition, dolomite was observed with XRD in the deepest sample of the northern-most pedon but not in pedons 13 or 22 (i.e. in those latter pedons, dolomite < 2 or 3 wt %).Williams (2008) developed a mineral model to calculate mineral abundances from bulk elemental analyses for the samples from pedons 1, 13, and 22 that matched the mineralogical data derived from quantitative X-ray diffraction.The mineral model included montmorillonite, illite, plagioclase, potassium feldspar, quartz, dolomite, and kaolinite plus minor calcite, chlorite, ferrihydrite.Using this mineral model, we calculated the mineralogy of the bulk elemental data from Muhs et al. (2001) for all the pedons along the north-south transect.Based on these calculations, the extent of depletion of plagioclase feldspar along the transect was observed to increase from north to south (Williams et al., 2010).In addition, evidence for a dolomite reaction front at depths between 200 and 400 cm was strong for pedons 1 through 6 (Mammoth, IL) and 8 (Greenbay Hollow, IL).The model was also consistent with a small amount of dolomite at depth in pedons 9, 17, 20-22.Consistent with this, other researchers have observed dolomite in deep loess samples near pedons 21 and 22 (Ruhe, 1984;Pye and Johnson, 1988;Muhs et al., 2001).To ascertain definitively if dolomite was present in any of these more southern profiles, we attempted to get samples in collaboration with D. Muhs (US Geological Survey).Sample volume limitations precluded analysis of most of the deepest samples from the southern part of the transect.In agreement with the mineral model, however, the one sam-ple (from 350 cm depth in pedon 20) that was large enough for micro-XRD revealed the presence of dolomite.
Given these observations, here we chose to follow Goddéris et al. (2010) by modeling pedons 1 in the north and 20 in the south.In both cases, we used mineralogy based on the mineral model of Williams (2008) (Table 1).This model predicts a small concentration of calcite; however, Godderis et al. (2010) documented that calcite removal in the WITCH model was rapid.Given that such small amounts of calcite might simply reflect error in the model attribution of mineralogy, we removed all calcite from the starting pedon mineralogy.For both pedons, the initial porosity was set equal to 43 % (Bettis et al., 2003).Furthermore, because the samples reported by Muhs et al. (2001) only allowed us to estimate that the dolomite reaction front was located between about 200 cm and 350 cm depth in pedon 20, whereas it was located at approximately 350 cm in pedon 1, we decided to fix the dolomite weathering front at 280 cm depth for both pedons.By locating the two fronts at the same depth, we can compare weathering in the north and south directly.Furthermore, moving this front by a few tens of cm up or down is not critical for this study.The front was also assumed to be a sharp front, meaning that we assume the total absence of dolomite above this depth, and an abundance of 11 % vol below this level (i.e.dolomite comprises a completely developed reaction front; Brantley and White, 2009).

Complex behavior of the driving parameters of weathering
It is generally assumed that weathering increases with runoff (runoff being the difference between rainfall and evapotranspiration) and with air temperature, at least for silicate rocks (White and Blum, 1995;Oliva et al., 2003;Dessert et al., 2003).Our simulation shows that the link between climate, vegetation and weathering is complex, especially when considered over a short timescale (here centuries).For example, the overall temperature rise from 4 to 6 • C calculated for the 1950-2100 period along the north-south transect, as predicted by ARPEGE, should promote silicate dissolution and slow carbonate dissolution because carbonates are generally at equilibrium with soil solutions and, unlike silicates, their solubility decreases with temperature (Schott et al., 2009).
However, increases in temperature are accompanied by a change in rainfall which is likely to also play a role.For the simulations, rainfall was predicted to rise moderately from 120 to 140 mm month −1 in the southern part of the Mississippi valley (31 • N) over the same time period, and from 75 to 85 mm month −1 for the northern area (42 since it also depends on changes in the actual evapotranspiration, which responds in a complex fashion owing to the rise in temperature and atmospheric CO 2 (including fertilization effects, stomatal closure, etc.).In the south, the yearly averaged vertical drainage calculated at the base of the pedon increases by 14 % (from about 350 mm yr −1 ) to a value slightly above 400 mm yr −1 -over the simulation period .In contrast, the deep drainage decreases in the north by more than 40 % (from 170 to < 100 mm yr −1 ) despite an increase in rainfall.Drainage decreases in the north because evapotranspiration (AET) increases significantly as temperature increases by 6 • C (Fig. 1).
In addition to changes in temperature and precipitation, the soil CO 2 concentrations also change in a complex fashion as global temperature and atmospheric CO 2 increase.As calculated by the CARAIB model, the temperature rise promotes soil respiration, but the higher temperature also accelerates CO 2 diffusive loss to the atmosphere (Van Bavel, 1951).Therefore, even though the modelcalculated soil CO 2 production rate rises from a mean annual value of 95 g C m −2 yr −1 to > 130 g C m −2 yr −1 in the south, the soil CO 2 concentration increases only from ∼ 20 × to 22 × atmospheric levels.This lack of change in the south is the outcome of faster CO 2 diffusion to the atmosphere that competes with the CO 2 fertilization effect, which promotes net primary productivity and hence both autotrophic and heterotrophic respiration.In contrast, this compensation does not occur in the northern site where the temperature effect on diffusion outweighs the fertilization effect causing soil CO 2 concentrations to decline from ∼ 30 × to 20 × atmospheric levels over the 1950-2100 period.This CO 2 decline happens despite a calculated increase in the soil CO 2 production rate from 70 to 110 g C m −2 yr −1 (Fig. 2).

CO 2 consumption by dolomite weathering
The instantaneous atmospheric CO 2 consumption by weathering is calculated as the flux of HCO − 3 required to precisely balance the cation release by the dissolution of silicate minerals, integrated over the whole weathering profile.In addition, the occasional precipitation of smectite contributes negatively to this budget.Finally, the bicarbonate flux equal to the Ca 2+ + Mg 2+ release by dolomite dissolution must also be added to account for atmospheric CO 2 consumed by dolomite weathering.
Despite the increase in temperature, the total CO 2 consumption due to weathering decreases in both the north and the south from 1950 to 2100.In the south, the model predicts a slight decrease from 1 to 0.9 mol m −2 yr −1 (trend calcu-  lated from a root mean square linear fit of the yearly averaged modeled CO 2 consumption).In the north, CO 2 consumption goes down by a factor of 2.3, from 0.7 to 0.3 mol m −2 yr −1 (Fig. 3).Dolomite dissolution is the main weathering reaction that largely controls these net CO 2 consumption trends: the CO 2 consumption flux by dolomite dissolution is about 10 × greater than consumption by any silicate phase in the pedons.In all the pedons, dolomite is absent at shallow depths and present only beneath 2.80 m depth.This 2.80 m depth therefore represents the sharp reaction front where downward-flowing fluids that are undersaturated with respect to dolomite begin to dissolve this mineral.The sharp reaction front and the close-to-equilibrium porefluids with dolomite are consistent with transport-controlled dissolution.Changes in the integrated dolomite weathering rate (in moles of mineral dissolved per time unit, summed over the whole weathering profile) are strongly correlated to the vertical drainage at both locations.In fact, as drainage increases over any time interval, dolomite dissolution responds linearly.This behavior is explained by the fact that waters draining the pedons are equilibrated with respect to dolomite below the dolomite front.Thus, the dissolution flux of dolomite is essentially equal to the product of the solubility of dolomite and the drainage flux (Brantley and White, 2009).If drainage increases, the flux increases, i.e., more water implies more dissolution.On the other hand, yearly averaged dolomite dissolution is negatively correlated to the temperature as dolomite solubility decreases with temperature.Since temperature displays a clear increasing trend from 1950 to 2100 period (4 • C warming) while drainage only increases slightly in the south or decreases drastically in the north (Fig. 1), the overall effect is a decrease in the CO 2 consumed by dolomite dissolution at both locations.

CO 2 consumption by primary silicate phases
In contrast to dolomite, the CO 2 consumption by dissolution of the primary mineral phases albite and K-feldspar either increase or stay at roughly constant values from 1950 to 2100.
In the south, equilibrium with respect to albite in the soil solution is episodically reached below 1 m depth, and permanently reached below the dolomite front at ∼ 2.80 m depth.
In the plots of soil solution saturation state with respect to albite versus time, pulses of high drainage are observed to be interspersed with periods of low drainage.High-drainage spikes lead to dilution below equilibrium in pore fluids down to the dolomite front, but the soil solution remains saturated with respect to albite below 2.80 m depth (Fig. 4a).The rate of occurrence of the dry events (low drainage periods when albite may effectively stop dissolving high in the profiles) increases after 2000 (Fig. 4a).But at the same time, the annually averaged drainage slightly increases in the south over the 1950-2100 period (14 % increase), promoting albite dissolution.The temperature rise (4 • C) forces albite dissolution to further increase.This results in an increasing trend in the calculated Na + export for the pedon from about 380 mol ha −1 yr −1 (1950) to 470 mol ha −1 yr −1 (2100).
For the northern pedon, the soil solution gets saturated with respect to albite mainly at the dolomite front even before 2000, i.e., at 2.80 m depth (Fig. 4b).This saturation level is reached at shallower depth after 2000.Several drought events strongly inhibit the vertical drainage, and these events become more and more common after 2000, owing to the overall decrease in drainage.As a result, the models predict an almost permanent rise of the saturation depth by about 2 m after 2030.However, this line cannot rise above the rooting zone because the through-flux of water is higher in that zone, as water that enters at the land surface is taken up by roots increasingly with depth.In contrast, beneath the rooting zone, the water flux is smaller and the water quickly equilibrates with respect to the low-solubility feldspar minerals (Fig. 4b).Like the dolomite behavior described above, the albite dissolution flux varies directly with vertical drainage in the north and in the south.In other words, just as dolomite dissolution is driven by the drainage flux, because dolomite is maintained at chemical equilibrium at the base of the pedons, the deep-soil solutions are at chemical equilibrium or even oversaturated with respect to albite.Any increase in drainage thus enhances the observed dissolution flux.In the north, the calculated decrease in drainage over the next century compensates for the temperature rise and the yearly averaged albite dissolution fluxes stay constant (or even slightly decrease) from 1950 to 2100.In contrast, in the south, where the drainage flux increases over the simulation period, the albite dissolution flux rises as the global temperature rises.
Similar arguments can be made for the behavior of Kfeldspar dissolution.For both pedons, soil solutions are saturated with respect to K-feldspar below the root zone at about 1 m.In the northern profile, dissolution increases from 1950 to 2100 and CO 2 consumption rises by 15 %, whereas in the south it rises by 25 %.
Not surprisingly, the flux from K-feldspar dissolution is positively correlated with vertical drainage at both sites.Indeed, like many two-feldspar systems (White et al., 2001), the saturation states of the draining waters with respect to K-feldspar are much higher than the saturation state with respect to albite.As a consequence, the dissolution flux from K-feldspar depends directly on the availability of water because everywhere in the profile dissolution occurs close to equilibrium with respect to the mineral phase: K-feldspar shows transport-controlled behavior.

CO 2 consumption by secondary silicate phases
Draining waters of both northern and southern sites are generally undersaturated with respect to smectite.Equilibrium with this phase is only attained during dry events in the northern site.For this site, the saturation state of the draining waters generally increases from 1950 to 2100 because of the occurrence of dry events related to the decreasing rainfall (Fig. 5).These dry events have a strong impact on the overall CO 2 consumption.As smectite precipitates, base cations are removed from the solution, limiting the atmospheric CO 2 consumption.Furthermore, these dry events reduce dolomite dissolution because drainage decreases down to 0 mm yr −1 below the root zone (Fig. 1b).As a consequence, CO 2 consumption falls to 0 mol m −2 yr −1 during extremely dry events such as the one indicated around 2025 (Fig. 3b).
The aridity developing in the north of the Mississippi valley promotes a general increase in the pH of the soil solutions (Fig. 6) (from 6.5 to 6.8 over 150 yr).There is a strong positive correlation between temperature and mean pH of the weathering profile.The correlation between mean annual drainage and mean annual pH is conversely negative, i.e., low-draining periods are characterized by higher pH.Both the temperature rise and the drainage decrease from 1950 to 2100 thus promote a pH rise in the soil solutions in the north.On the contrary, pH oscillates around an average value of 6.65 over the whole simulation period for the south pedon.It is also interesting to note that porewaters available to plants change in composition and pH forward in time.In other words, the chemistry at rooting depth is affected by climate change and plants have access to fluids of different chemistry over time.Of course, there is no feedback of the soil chemistry on the biospheric model in the present study, and such feedback should be considered in the future.

Relative contribution of the dissolution of each mineral to CO 2 consumption
Overall, the relative contribution of each mineral phase to total CO 2 consumption depends largely on the vertical drainage at both sites.The dolomite contribution to total CO 2 consumption is generally above 90 % but can decrease to less than 40 % during dry events in the north (Fig. 7).
Conversely, the relative silicate contribution increases during low drainage events.Indeed, the associated CO 2 consumption flux goes down when drainage decreases, but at a slower rate than the flux due to dolomite alone.

Calculated CO 2 consumption compared to previous estimates
Our simulations predict a CO 2 consumption of about 1 mol CO 2 m −2 yr −1 for the southern pedon and 0.3-0.7 mol CO 2 m −2 yr −1 for the northern pedon.Since most of the CO 2 consumption is the result of dolomite dissolution, these numbers can be compared to the CO 2 consumption by outcropping carbonate globally.Carbonate outcrops (including marls, dolomites, and limestones predominantly) cover about 13.4 % of the continental surface (20.1 × 10 6 km 2 ) (Amiotte- Suchet et al., 2003).They consume 12.3 × 10 12 mol of atmospheric CO 2 yr −1 (Gaillardet et al., 1999).This gives a mean area-averaged CO 2 consumption of 0.62 mol m −2 yr −1 , close to our estimated values for the loess.
Results can be also compared with the modeled values (using parametric laws calibrated on field data) from Moosdorf et al. (2011).For unconsolidated sediments of North America, including loess, CO 2 consumption is estimated at 0.2 mol CO 2 m −2 yr −1 for the southern pedon using the Moosdorf et al. (2011) correlation between CO 2 consumption and runoff (taken as the difference between rainfall and evapotranspiration), and about 0.05 mol CO 2 m −2 yr −1 for the northern pedon.For carbonate rocks of North America, these numbers rise to 0.8 and 0.25 mol CO 2 m −2 yr −1 , respectively, in Moosdorf et al. (2011).These numbers are of the same order of magnitude as our results.

Sensitivity tests
The key role of temperature on dolomite dissolution and the associated trends in CO 2 consumption can be illustrated by a simulation where monthly air temperature is fixed at its monthly averaged 1950-1960 value.In the south, the overall CO 2 consumption now increases at a secular rate of 1 mmol CO 2 m −2 yr −2 , while it was decreasing at a rate of −1 mmol CO 2 m −2 yr −2 in the reference run.In the north, the decreasing trend of −2.5 mmol CO 2 m −2 yr −2 (reference run) is converted into an increasing rate of 0.25 mmol CO 2 m −2 yr −2 (sensitivity test; Fig. 3).
This does not mean that temperature is a stronger driver of the CO 2 consumption by weathering than drainage for these soil profiles.Indeed, for both sites, the base cation export fluxes at the base of each pedon are strongly correlated to the vertical drainage.But the vertical drainage only slightly increases in the south while it strongly decreases in the north, owing to enhanced evapotranspiration in a warmer world.The response in terms of weathering is strongly dependent on the climate and biospheric model behavior: temperature appears as a key factor because the ARPEGE GCM predicts a rather strong global warming in the Mississippi valley (4 to 6 • C), which impacts the biospheric model in such a way that the evapotranspiration increase compensates for the rainfall increase, limiting the drainage response.

Carbonate fronts along the Mississippi transect
The detailed analysis presented above was focused on the most southern and northern pedons because those are the two sites where a precise mineralogical composition has been calculated as a function of depth.Nevertheless, these two mineralogical profiles are not widely different (see Table 1).We thus expand our calculations to the whole Mississippi latitudinal transect, assuming the mineralogical composition of the southern pedon for all sites.As discussed by Williams et al. (2010), the Peoria loess is fairly invariant in composition.Our goal is to explore the response to climate change of the most reactive mineral phase with respect to the CO 2 budget (here dolomite, as demonstrated above), along a latitudinal transect.
The dolomite front corresponds to the depth below which soil solutions are supersaturated with respect to dolomite, and above which they are undersaturated.The dolomite reaction front can be viewed as a "terrestrial lysocline", i.e. a depth in the soil where the rate of dissolution of carbonate changes rapidly, like the oceanic lysocline.But contrary to the downward-flowing soil solutions, seawater is oversaturated with respect to carbonate phases (either calcite or aragonite) above the oceanic lysocline and undersaturated at depth.With the ongoing oceanic acidification process driven by the rising CO 2 , the oceanic lysocline is moving upward, reducing the thickness of the ocean layer saturated with respect to carbonate minerals (Archer et al., 1997;Boudreau et al., 2010).
However, the two lysoclines are coupled.When the ocean lysocline shallows, the terrestrial lysocline is deepening such that after some time lag, a new steady-state ocean pH and atmospheric CO 2 concentration must be achieved.In effect, the enhanced continental carbonate dissolution produces a flux of alkalinity to the oceans which eventually can counteract acidification (Archer, 2005).In other words, under some conditions, the depth of carbonate depletion on the continents would increase to compensate for the thinning of the carbonate stability zone in the oceans.Such a feedback is common in studies of the geological evolution of the climate, where the surficial ocean stays oversaturated with respect to carbonate despite atmospheric CO 2 levels reaching 10 times or more the present-day level.The question is whether this feedback can work at the secular or millennial timescale.The answer to this question is not obvious, since carbonate dissolution (here dolomite) is heavily dependent on the temperature and drainage evolution into the future, and this latest parameter displays contrasting behavior along the Mississippi transect.
Model calculations for pedons over the entire transect document that the dolomite front deepens 1.9 cm over the 1950-2100 period in the northern part of the Mississippi (40.9 • N), and 6.1 cm over the same time period in the central part of the transect (34.2 • N) (the greatest value observed for the transect).These rates correspond to a retreat rate of 13 cm kyr −1 in the north, and 41 cm kyr −1 in the central part.The amplitude of the retreat is approximately constant (≈ 5 to 6 cm retreat over 150 yr) from 30 • N to 37 • N. Above 37 • N, the dolomite front retreat rate suddenly decreases below 2 cm over the 150 yr of the study.This sudden drop correlates well with a marked decrease in rainfall around 37 • N. The averaged rainfall over the whole period of simulation falls from about 105-130 mm month −1 below 37 • N to 80 mm month −1 above this latitude.
This rate of retreat of the front is correlated to the local climatic conditions, but does not carry any information about the response of the terrestrial lysocline to climate change.But the question remains: does the rate of retreat accelerate or not over the next century?We first calculated the time derivative of the depth of the terrestrial lysocline for 9 pedons.These rates of retreat display ample fluctuations as a function of time for each location.We isolate the long-term trends by fitting these rates (in cm yr −1 ) with a linear fit.The slope of this linear fit (cm yr −2 ) is a measurement of the long-term acceleration or deceleration of the retreat rate of the terrestrial lysocline (Fig. 8).Over the 1950-2100 period, the acceleration of the retreat rate is maximum in the central part of the transect.In contrast, it is at a minimum, reaching a negative value, in the north, suggesting that the dolomite front is still retreating downward, but at a rate that decreases over the 1950-2100 period.In the south, the downward shift of the dolomite front is decelerating, but at a moderate rate.
The causes for this contrasting behavior are multiple.In the north, temperature is increasing at the fastest rate for any of the sites along the transect, soil CO 2 is decreasing rapidly, and drainage is decreasing.All these factors tend to slow down the dolomite dissolution and hence the downward rate of retreat of the terrestrial lysocline.In the central part of the watershed, in contrast, the drainage increase is at a maximum, accelerating the rate of retreat of the dolomite front.In the south, processes are more intertwined, with the rate of increase in soil CO 2 playing an important role for this environment where the changes in drainage and temperature are similar.It is important to note that the drainage rise in the extreme south is not large enough to compensate for the temperature rise.As a consequence, the dolomite front is retreating there at a decreasing rate for the 1950-2100 period.
We note again that the initial position of the dolomite front was fixed at 2.8 m depth for all pedons.Moving this initial position upward or downward does not alter significantly our calculated retreat rates.Indeed, dolomite front retreat rates are mostly affected by changes in pH (Fig. 6) and drainage, and these changes largely occur within the root zone.pH and drainage remain fairly constant with depth below the root zone (until the dolomite front where pH rapidly rises).As long as the dolomite front is located below the root zone, which is the case for all pedons according to Williams et al. (2010), the geographical pattern of the terrestrial lysocline retreat rate is not dependent on the initial dolomite front depth.
These factors (drainage, temperature and soil CO 2 ) thus produce complex behavior that cannot easily be integrated or averaged over the transect to answer the question we have asked.In the central part of the watershed, dolomite weathering is accelerating, while it is slowing rapidly in the north and moderately in the extreme south.The response of the terrestrial lysocline to the ongoing climate changes is not straightforward.Indeed, Fig. 8 illustrates the high sensitivity of the terrestrial lysocline deepening to parameters which are themselves sensitive to the modality of climate change (drainage, temperature and soil CO 2 levels).Fig. 8. Increasing rates of vertical drainage (below the root zone), ground temperature, and soil CO 2 level over the Mississippi transect, simulated for the 1950-2100 period.In red, acceleration of the dolomite retreat rate along the transect simulated over the same period.The dolomite front is moving downwards at each location, but this retreat rate either accelerates (positive value) or decelerates (negative value).

Y. Goddéris et al.: Rates of consumption of atmospheric CO
Maximum acceleration of the retreat rate is reached in the central section of the Mississippi Valley, owing to the rate of change in climatic and geochemical conditions.

Limitations and perspectives
Contrary to the temperature and soil hydrology, the soil CO 2 level does not emerge amongst the critical driving factors of weathering in the future in our simulations.However, our model for soil CO 2 is rather basic.Diffusion only depends on air temperature (Van Bavel, 1951;Goddéris et al., 2010).The impact of soil hydrology on diffusion is not included, e.g. the expected decrease in air-filled porosity with increasing water content, which would limit diffusion back to the atmosphere when the profiles get more humid.A better model should be implemented in the future.However, a striking output of the simulation is the decrease in soil CO 2 in the future for the northern pedon, as diffusion rises with temperature.As the simulated profile also dries out, upward diffusion would be further promoted.Thus the soil CO 2 trend will still display a decrease into the future even if the role of hydric state on diffusion is included.Another limitation arises when considering the terrestrial lysocline behavior.In our simulation, we do not allow dolomite to precipitate when solutions become supersaturated.Indeed, the kinetics of dolomite precipitation is very low at ambient temperature.Globally, dolomite precipitation is usually only observed where promoted by bacterial activity in specific environments (Vasconcelos et al., 1995;Sánchez-Román et al., 2011), but calcite is not allowed to precipitate either.Supersaturation might occur in the future, especially since the weathering profiles get warmer and possibly dry out, as shown in the simulations.If calcite were to precipitate, the terrestrial lysocline might cease deepening and actually move upward, consuming alkalinity and releasing CO 2 back to the atmosphere.Although atmospheric carbon would be trapped in mineral form if calcite precipitates, there would be less capture of carbon in the dissolved HCO − 3 form.This latter capture stores carbon in the world oceans at the millen-nial timescale (Beaulieu et al., 2012).An upward movement of the terrestrial lysocline would negatively impact the short term atmospheric CO 2 consumption by weathering.
Here, below the dolomite front in both north and south sites, the saturation state of the soil solution with respect to calcite is predicted to increase in the future (Fig. 9).In the south, the mean annual saturation state of the soil solution with respect to calcite below the dolomite front is seldom above 1 before about 2015.Then the occurrence of saturation increases rapidly after 2015.In contrast in the north, saturation of the soil solution stays generally below 1, because temperature is lower, although it increases with time.
Although not critical here, future work should account for calcite precipitation.
Finally, our work shows that the weathering of carbonate lithologies may respond quickly to climate change in the near future.Because of the very high reactivity of carbonates, their dissolution is mainly controlled by transport, stressing the need for an accurate modeling of the water fluxes in these environments.Modeling of carbonate precipitation, which may occur more often in a drier and warmer future world, would require accurate knowledge of the carbonate precipitation kinetics close to equilibrium.

Conclusions
We use a cascade of numerical models (simulating climate, the continental biosphere, and the weathering processes) to extrapolate the weathering of the Mississippi Valley loess forward in time to 2100.Our interest in the weathering of those loessic formations arises because they are globally important due to their high surface area and due to the fact that they contain both carbonate and silicate mineralogies, thus exemplifying the contrasting behaviors of weathering  for these two dominant rock-forming mineral types.However, loess is also of importance in that it represents a fastweathering component of the terrestrial surface that is produced in response to glaciations worldwide (Anderson et al., 1997).The simulations we present here are an outgrowth of our previous investigation of the Peoria loess deposited along the Mississippi valley of the USA interpreted as a natural experiment for pedogenesis over the last 13 kyr (Williams et al., 2010;Goddéris et al., 2010).
We found that yearly averaged drainage is expected to slightly increase in the southern end of the Mississippi Valley (14 % increase), while it will decrease drastically in the northern end (40 % decrease).The decrease in the north is triggered by the rise in evapotranspiration driven by enhanced temperatures.Coupled to the rise in temperature linked to the rise in CO 2 (700 ppmv in 2100), the dolomite dissolution rate slows down in the northern and southern ends of the Mississippi Valley.Dolomite dissolution is further reduced by the accelerated diffusion of soil CO 2 to the atmosphere under a warmer climate.Conversely, primary feldspar dissolution rates are enhanced by the temperature rise.
Dry events occur more often after about 2025, particularly in the north, forcing pH to rise by 0.3 units in the weathering profile.Smectite precipitates during very low drainage events.
When simulating the whole Mississippi transect (9 sites, from 30.7 • to 42.2 • N), we found that the dolomite front retreat decelerates at both ends of the transect, but accelerates in the middle part of the watershed.This area experiences the strongest enhancement in vertical drainage as dictated by the delicate balance between rainfall and evapotranspiration along the transect.In the middle part of the watershed, this balance maximizes the drainage increases.North of 37.4 • N, the drainage rapidly decreases, forcing the dolomite retreat rate to slow down into the future.
The dolomite front defines a terrestrial lysocline for the Mississippi Valley loess.The deepening of this lysocline defines the way carbonate weathering responds to climate change.The deepening accelerates only in the central part of the watershed, and decelerates elsewhere.Acceleration of the retreat of the terrestrial lysocline is equivalent to increasing rates of delivery of alkalinity to the oceans.Defining how this lysocline will respond at continental or global scales is thus a key issue regarding predictions of impacts of processes such as ocean acidification that are occurring as a result of increasing atmospheric CO 2 .

Model details A1 WITCH general description
WITCH is a vertical box model that calculates the time evolution of the soil solutions as a function of time.Soil solution chemical composition is calculated for 40 vertical levels extending from the surface down to 4 m depth.The vertical resolution is 0.2 m in the upper horizons, and decreases below the dolomite weathering front to match the fluctuations in the advance rate of this front correctly.A differential equation describes the budget of the main conservative species (not involved in fast speciation reactions) and is solved at each time step for each box: where C j is the concentration of the species j .WITCH includes a budget equation for Ca 2+ , Mg 2+ , K + , Na + , SO 2− 4 , total alkalinity, total aluminum, total silica, and total phosphorus.z is the thickness of the layer under consideration, Y. Goddéris et al.: Rates of consumption of atmospheric CO 2 and θ is the water volumetric content.F top is the input of species from the layer above the given layer in moles per unit time through drainage, and F bot is the removal of water through downward drainage.n m is the total number of minerals considered in the simulation.R j is the exchange term between the soil solution and the clay-humic complex.F i weath, j is the release of species j to solution through dissolution of the mineral i in the box, and F i prec, j is the removal of j through precipitation when the solution becomes supersaturated with respect to mineral i.The dissolution F i weath, j and precipitation F i prec, j terms for silicate minerals are calculated using kinetic laws and parameters derived from transition state theory (TST; Eyring, 1935) and laboratory experiments, respectively (Schott et al., 2009;Goddéris et al., 2010): where F g is the overall dissolution rate of mineral g inside a given layer.The sum accounts for the four parallel ratecontrolling elementary reactions that are assumed to describe dissolution promoted by H + , OH − , H 2 O and organic ligands.
In these relations, a l and n l, g stand for the activity and the reaction order with respect to the i-th species, respectively.For organic ligands, a l equals the activity of a generic organic species RCOO − (conjugate of RCOOH).k l, g is the rate constant of mineral g dissolution reaction promoted by species l, and E l a, g is the activation energy of this reaction.f inh stands for inhibitory effects (i.e. by aqueous Al).The last factor in Eq. (A2) describes the effect on the rate of the departure from equilibrium, where g is the solution saturation state with respect to mineral g ( g = Q/K g where Q is the activity quotient and K g is the equilibrium constant for mineral g dissociation reaction), and s, the Temkin's coefficient number, is the stoichiometric number of moles of activated complex formed from one mole of the mineral.Precipitation terms are calculated in the same way, but g is then > 1 (only clay secondary phases are allowed to precipitate).A g is the reactive surface.Dolomite dissolution has been added to WITCH and modeled within the framework of TST and surface coordination chemistry concepts (Pokrovsky andSchott, 1999, 2001;Pokrovsky et al., 1999).Dolomite dissolution rate F dol is modeled as: 1.575 • 10 −9 1.575 • 10 −9 + 3.5 • 10 −5 • a CO 3 + a CO 3 • a Ca • 1 − 1.9 dol ,
Once the mass balance is calculated, the speciation of the soil solutions is calculated at each time step together with pH.In addition, WITCH includes a mass balance for each mineral phase, accounting for the amount of mineral being dissolved and precipitated within each time step.

A2 Soil hydrology
The water volumetric content θ at each weathering profile level and the vertical drainage as a function of depth are taken for each location from the simulation of the global dynamic vegetation model CARAIB, as mentioned in the main text (Dury et al., 2011).CARAIB accounts for only one level below ground, spanning the whole root zone, while the vertical resolution of WITCH is much finer.Regarding θ, we fix the water volumetric content of the first layer of WITCH below ground at the water volumetric content calculated by CARAIB.At the base of the root zone (a depth which is calculated at each time step by CARAIB), we fix θ at the field capacity.Inside the root zone, the water volumetric content is linearly interpolated.
Regarding the vertical drainage, the water flux at the ground level and entering the weathering profile is fixed at the rainfall calculated by the ARPEGE GCM minus the direct evaporation.At the base of the root depth, the drainage in WITCH is fixed at the drainage calculated by CARAIB, accounting for the water uptake by the roots.Below that level, drainage is held constant.Inside the root zone, drainage is linearly interpolated at each time step.Apart from the surficial runoff at the top of the weathering profile, the model does not account for horizontal transfer inside the profile.Once entering the profile, the water is either taken up by roots or transferred downward.

A3 Calculating soil CO 2 levels
The CARAIB model calculates autotrophic and heterotrophic respiration.We assume that all the heterotrophic respiration and 1/3 of the autotrophic respiration occurs below ground.The sum of these two terms defines the CO 2 production inside the root zone, the thickness of which is calculated by CARAIB.This production allows calculation of soil CO 2 at the base of the root zone, accounting for the temperature effect on CO 2 upward diffusion, as detailed in Goddéris et al. (2010).We then assume that the CO 2 pressure at ground level (0 m depth) equals the atmospheric pressure.CO 2 is assumed to increase linearly inside the root zone until it reaches its maximum at the base of the root zone.Below this level, it is held constant.
Edited by: J. Middelburg , among Published by Copernicus Publications on behalf of the European Geosciences Union.Y. Goddéris et al.: Rates of consumption of atmospheric CO 2

Fig. 1 .
Fig. 1.The calculated water budget for the southern ((a) and (b)) and northern ((c) and (d)) sites.The increase in evapotranspiration driven by the temperature rise offsets the drainage response to the global warming.The 1-yr running means of the monthly outputs are plotted.The dotted lines are the mean square linear fits of the vertical drainage evolution at the basis of the pedons.

Fig. 4 .
Figure 4 Fig. 4. Calculated time evolution of the saturation state of the soil solutions with respect to albite, plotted as a function of depth: (a) southern pedon, (b) northern pedon.

Figure 5 Fig. 5 .
Figure 5 Fig. 5. Calculated smectite dissolution rates plotted as a function of depth and time for the northern pedon.

Fig. 6 .
Fig. 6.Calculated soil solution pH plotted as a function of depth and time for the northern pedon.

Figure 7 Fig. 7 .
Figure 7 Fig. 7. Relative contribution of calculated dolomite dissolution to the CO 2 consumption in the north plotted as a function of vertical drainage.Each point stands for a yearly averaged value.

Fig. 9 .
Fig. 9. Calculated saturation state of the soil solutions with respect to calcite plotted as a function of depth and time: (a) south pedon, (b) north pedon.
k dol H equals 10 −3 mol m −2 s −1 , and k dol Mg 10 −8.2 mol m −2 s −1 at 25 • C. The activation energy for k dol H and k dol Mg is set to 40 and 60 kJ mol −1 , respectively (Pokrovsky

Table 1 .
• N).The vertical drainage through the loess formation will also vary as rainfall increases.However, drainage does not vary in a simple way Pedon mineralogy (%) assumed for pedon north (42 • 59 N) and pedon south (30 • 47 N).