Global weathering of calcium and magnesium silicate rocks provides the long-term sink for atmospheric carbon dioxide (CO2) on a timescale of millions of years by causing precipitation of calcium carbonates on the seafloor. Catchment-scale field studies consistently indicate that vegetation increases silicate rock weathering, but incorporating the effects of trees and fungal symbionts into geochemical carbon cycle models has relied upon simple empirical scaling functions. Here, we describe the development and application of a process-based approach to deriving quantitative estimates of weathering by plant roots, associated symbiotic mycorrhizal fungi and climate. Our approach accounts for the influence of terrestrial primary productivity via nutrient uptake on soil chemistry and mineral weathering, driven by simulations using a dynamic global vegetation model coupled to an ocean–atmosphere general circulation model of the Earth's climate. The strategy is successfully validated against observations of weathering in watersheds around the world, indicating that it may have some utility when extrapolated into the past. When applied to a suite of six global simulations from 215 to 50 Ma, we find significantly larger effects over the past 220 Myr relative to the present day. Vegetation and mycorrhizal fungi enhanced climate-driven weathering by a factor of up to 2. Overall, we demonstrate a more realistic process-based treatment of plant fungal–geosphere interactions at the global scale, which constitutes a first step towards developing ‘next-generation’ geochemical models.
On geological timescales, atmospheric carbon dioxide (CO2) is largely controlled by the following generalized reaction [1–3]: 1.1which describes the weathering of calcium-bearing silicate rocks on land by atmospheric CO2 followed by the deposition of carbonate minerals on the ocean floor. Release of magnesium from silicate rocks during weathering also leads to limestone deposition following the exchange of magnesium for calcium on contact with ocean-floor basalts , whereas other base cations (K, Na) do not contribute to the large-scale deposition of carbonates. Silicate minerals contain no CO2, so release of their calcium and magnesium cannot lead to precipitation of limestone without removing CO2 from the atmosphere–ocean–soil system. Because the limestones will remain stable on the ocean floor until tectonic processes exhume or destroy them, equation (1.1) depicts the sequestration of CO2 on timescales of millions of years. Equation (1.1) forms the basis of geochemical numerical models of the long-term carbon cycle [3,5–10], which calculate atmospheric CO2 as the balance between weathering and degassing on multi-million-year timescales [4,11].
Vegetation has long been thought to have exerted a profound influence on climate through weathering [12–14] by altering soil hydrology and stability, and the pH and organic ligand content of the soil solution. This idea has recently been extended to recognize the influence of mycorrhizal fungi , which profoundly alter the nature of plant–mineral interactions, directing the host plant's photosynthetic weathering energy to the surfaces of nutrient-bearing minerals more effectively than roots. The ancestral arbuscular-mycorrhizal (AM) fungi (Glomeromycota) have behaved as biosensors, seeking nutrients and increasing the effective surface area of contact between minerals and the root systems of plants, since at least Devonian times . The more recently evolved ectomycorrhizal (EM) fungi (Basidiomycota and Ascomycota) further enhance this effect by effectively preventing root tip contact with soil, so that the majority of nutrient uptake is undertaken by the fungus.
Building on this conceptual framework, Taylor et al.  developed a process-based model investigating the possible role of plants and mycorrhizal fungi in biological weathering of silicate minerals. The model describes how the primary productivity of terrestrial vegetation and associated nutrient uptake by fine roots and mycorrhizal fungi influences soil chemistry and mineral weathering via its effects on the biological proton cycle, and integrates these processes into a leading global geochemical carbon cycle model (GEOCARBSULF). The biological proton cycle is controlled by the stoichiometry of reactions during the uptake of mineral nutrients by plants and symbiotic root fungi during growth, and the subsequent return of these nutrients during decomposition and leaching of organic matter . Overall, these reactions determine the net flux of protons into and out of the sub-surface environment and, therefore, the acid–base balance and its influence on the pH of the mineral weathering environment. The mechanistic analyses of Taylor et al.  indicated that the evolution and spread of EM fungi at the expense of the ancestral AM fungi could have contributed to CO2 drawdown in the Cretaceous. Although the authors assumed that the bulk of the silicate Ca + Mg weathering flux came from ‘hotspots’ , rather than being evenly distributed across the globe, they made no attempt to locate or calculate the influence of such hotspots on the global weathering flux.
Models of the long-term geochemical carbon cycle generally adopt a zero-dimensional approach, where the entire global land surface is reduced to one uniform box which weathers at the mean land surface global temperature. This is a useful first-order approach but cannot capture spatially explicit effects or the effect of changing continental configurations on temperature and runoff [18–21], nor the location of highly reactive volcanic rocks with respect to regional or local climate and vegetation. Therefore, we now extend the modelling approach of Taylor et al.  into two dimensions (latitude by longitude) to investigate the effects of regional and seasonal variations in the activity of vegetation and mycorrhizal fungi as well as climate on continental weathering processes over geological time. Our global strategy couples our biological proton cycle weathering model  with the Sheffield dynamic global vegetation model (SDGVM) , which simulates terrestrial carbon cycle processes, and the Hadley Centre coupled ocean–atmosphere general circulation model (HadCM3L) [23,24]. Simulated global patterns of vegetation activity, the geographical distribution of plant functional types (PFTs), land surface hydrology, climate and lithological maps are then used to simulate the biological proton cycle, soil solution chemistry and mineral weathering to calculate spatially resolved continental Ca + Mg fluxes.
Our approach parallels and significantly extends earlier attempts using dynamic global vegetation models (DGVMs) to estimate biological weathering, because we focus on the effects of biological activity on soil solution chemistry. Roelandt et al.  coupled a DGVM into a sophisticated n-layer mechanistic weathering model  requiring detailed information about soil characteristics for each grid point, that accounts for both evapotranspiration and root respiration effects on weathering via hydrology and soil pCO2. Unlike the simpler one-layer box model of Taylor et al. , it takes no account of the actions of mycorrhizal fungi on the biological proton cycle that control mineral weathering . It is important to note that this study focused on the catchment scale. Other investigations have attempted global-scale weathering calculations by coupling a DGVM with a model of the geochemical carbon cycle, but these focused on characterizing the indirect role of changing continental vegetation on weathering via its effects on climate by changing surface albedo, hydrology and land surface roughness .
We present this paper in three parts. First, we outline the development of our global-scale modelling strategy. Next, we report evaluation of the approach for contemporary ecosystems under modern climate and CO2 by comparing simulated weathering fluxes at the catchment scale with measurements worldwide . We use contemporary climate and simulated ecosystem properties to address the role of vegetation primary production, and plant and mycorrhizal functional type in driving continental weathering. Finally, we extend the modelling approach back in time with a suite of HadCM3L palaeoclimate simulations spanning the past 220 Myr to provide provisional estimates of continental-scale biological weathering over this time. Our aim here is to use a process-based approach to estimate how weathering by land plants and fungi may have changed over time relative to the present day. This is a critical response integrated into all geochemical carbon cycle models applied to the Phanerozoic [5,7,8,19,29], but is presently represented empirically in a manner that assumes the efficiency of weathering in the past is not greater than it is for modern contemporary ecosystems.
2. Model development
Figure 1 indicates that silicate mineral weathering regimes are influenced by a variety of biotic and abiotic factors. Plate tectonics affects regional climate via continental position and relief, and global climate through degassing of CO2, but the high relief, flood basalts and andesitic island arcs generated at plate boundaries also directly impact silicate weathering fluxes of Ca and Mg from land to oceans.
The three main components of our two-dimensional model (Hadley Centre CM3, SDGVM and our weathering model) are coupled as indicated in figure 1.
(a) Climate and vegetation models
HadCM3L is a new generation of coupled ocean–atmosphere models needing no flux correction terms that is well-documented and tested elsewhere [23,24]. This particular version uses the standard resolution atmosphere (3.75° × 2.5° × 19 levels) coupled to an equivalent resolution ocean (3.75° × 2.5° × 20 levels), with MOSES2.2 land surface scheme. It generates monthly average temperature, precipitation and humidity for each grid point, given external boundary conditions including orography, continental position and global atmospheric CO2. We forced HadCM3L with appropriate palaeogeographies provided by Markwick & Valdes . The solar constant was modified to be appropriate for the relevant period.
Terrestrial carbon cycle processes were modelled with SDGVM, which has been extensively described and validated against observations elsewhere . In brief, SDGVM simulates global patterns of net primary production (NPP), leaf area index (LAI) and the distribution of PFTs from monthly inputs of temperature, precipitation, relative humidity and global datasets of soil texture . Core modules of net photosynthesis, stomatal conductance, canopy transpiration, uptake of mineralized nitrogen and responses of these attributes to changes in soil water supply are detailed, and rigorously evaluated against field observations [31,32]. A key feature of SDGVM is the coupling of above- and below-ground C and N cycles. Litter production influences soil C and N pools via the Century soil nutrient cycling model , which feeds back to influence above-ground primary productivity. Local- and global-scale SDGVM predictions of NPP, LAI and PFT distributions have been extensively and successfully evaluated against a wide range of measurements, field observations and satellite products [22,32].
(b) Weathering model
As described by Taylor et al. , our weathering module assumes that soil pH in the mycorrhizosphere (i.e. the volume immediately around the root and mycorrhizal fungal hypha with an effective radius of several millimetres) is largely controlled by the ‘biological proton cycle’ describing the stoichiometry of proton and base cation cycling between soil, organisms and organic matter . Overall, these reactions determine the net flux of protons into and out of the sub-surface environment and therefore the acid–base balance and its influence on the pH of the mineral-weathering environment. Numerical routines account for weathering rates on basalt and granite that incorporate simple yet rigorous equilibrium chemistry  and rate laws [35,36]. Mineral dissolution reactions are described with kinetic mass action laws for irreversible reactions far from equilibrium. Aqueous speciation reactions such as inorganic carbon speciation and pH buffering are treated as thermodynamic equilibria. The model includes representation of the biological cycling of protons and alkalinity  occurring in discrete regolith zones: (i) the soil that is in contact with or under the influence of living roots and mycorrhizal hyphae (the mycorrhizosphere) and (ii) the remaining bulk soil. Weathering in both zones is influenced by abiotic processes and decomposition .
Our simulation of the influence of vegetation and fungi on the biological proton cycle at the global scale begins by considering the fine roots and hyphae required to take up sufficient nutrients to support the NPP simulated by SDGVM, assuming an average stoichiometry . Fine roots and EM hyphae also exude low molecular weight organic acids at rates of 10 mol m−1 fine roots day−1  and 10 mol m−1 hyphae day−1 , respectively. A further contribution of organic acids is linked to leached carbon provided by SDGVM, because removal of organic matter is associated with permanent loss of alkalinity from the soil [10,16]. Because SDGVM produces steady-state vegetation, we do not model the effects of aggrading ecosystems in this paper.
The soil depth distributions of fine roots and mycorrhizal hyphae are modelled using an exponential function , 2.1where ϕ is the length at depth z, ϕ0 is the length at the top of the soil and is a characteristic depth defined such that 2.2A similar expression describes the distribution of CO2 with depth, applying an analytical solution to the diffusion equation .
Fine roots and hyphae are assigned lengths appropriate for modern plant and mycorrhizal functional types [41,42]; these lengths and the associated characteristic depths are allowed to vary linearly as NPP deviates from modern averages for those same functional types (C3 and C4 grasses, and trees with any combination of deciduous or evergreen and broad or needle leaves). We can then calculate the integral of the distribution described by equation (2.2) to define the volume of the rhizosphere for the soil under 1 m2 of land. Assuming a uniform distribution of soil water, we can combine this soil volume with the soil water content provided by SDGVM to get the rhizosphere solution volume. All surfaces in the rhizosphere are influenced by nutrient uptake and organic anion exudation by fine roots. Mycorrhizal fungi, on the other hand, are assumed to seek the minerals containing the most available nutrients , and their uptake and exudation effects are restricted to those minerals. This means that hyphal interactions with minerals are amplified in comparison with fine roots because the hyphosphere is a sub-volume of the rhizosphere.
These descriptions of nutrient uptake, exudation and carbon loss, combined with the rhizosphere and hyphosphere solution volumes, estimates of the residence time of water in micropores and rate laws of the form 2.3allow us to calculate the alkalinity and thus the pH of the film of water on mineral surfaces . In equation (2.3), the subscripts i and j refer to minerals and weathering agents aj (H+, H2O, OH− and oxalate), respectively; Eapp,i is the apparent activation energy for weathering mineral i; kij and nij are the rate constant and reaction order for mineral i undergoing weathering by agent aj; R is the gas constant; and T is the temperature (K).
(c) Mineral surface area and lithology
We calculate weathering fluxes with rate laws for the silicate base-cation-bearing minerals found in basalt, granite, shale and sandstone, including feldspar and plagioclase, inosilicates, biotite and clay minerals . Only silicate minerals are included; carbonates are ignored because weathering of these minerals releases CO2 and, therefore they are less important to the long-term carbon cycle. The rate laws for H+, H2O and OH− are from the compendium of Palandri & Kharaka .
We approximate the fraction of the reactive surface area of interest for each mineral using the weight per cent of the mineral in the rock. The total reactive surface area under 1m2 of land is a fixed estimate for each rock type. For basalt, we use 20 000 m2 m−2 land, the scaling value derived by Navarre-Sitchler & Brantley  for comparing watershed-scale weathering rates to laboratory rates on basalt substrates. Our mineralogy is that of a normal tholeittic basalt  and includes the Ca- and Mg-bearing minerals labradorite, augite and orthopyroxene.
For granite, we use the mineralogy of the young soils in the chronosequence examined by White et al.  for which we derive a total surface area of 1 788 040 m2 m−3 soil for the silicate size fraction greater than 4 µm. Because changes in size fractions along their chronosequence indicate reaction rates four orders of magnitude lower than laboratory rates, we use a reactive silicate surface area of 179 m2 m−3 soil. This number also compares well with the estimate of 240 m2 m−3 soil for granite at 0 Ma by Taylor et al. . The minerals providing base cations are K-feldspar, plagioclase, hornblende and biotite. We calculated global average compositions for shales and sandstones from the platform and geosynclinal averages given by Holland , using the same procedure he used to derive an average global sedimentary composition. The fractions of Ca, Mg, CO2 and SO42− indicate that 20 per cent of Ca and 66.4 per cent of Mg in sandstones is found in silicate minerals; we assumed that all of the Na and silicate Ca is in albite (Ab92), all of the K is in feldspar, and the silicate Mg is in chlorite. Likewise, our shale has 28.1 per cent of its Ca and 71.9 per cent of its Mg in silicate minerals. Shales may include smectites with variable proportions of base cations; we used the composition for which Palandri & Kharaka  provide a rate law: K0.04Ca0.5(Al2.8Fe0.53Mg0.7)(Si7.65Al0.35)O20(OH)4. We assumed that all the Na was in the same albite we used for sandstone (Ab92), all the remaining silicate Ca was in smectite and all the remaining Mg in chlorite. We did not attempt to alter this mineralogy to reflect latitudinal variations in clay composition.
Brandt et al.  found that the reactive surface area of chlorite corresponded to molecular steps in the t–o–t layers comprising just 0.2 per cent of the surface area determined by gas adsorption. If a shale comprised entirely clay minerals with surface areas of 0.0024 m2 g−1 soil and specific gravity G = 2.6, then one can calculate a surface area of 6240 m2m−3 soil. Because clay minerals comprise only 24 per cent and 4 per cent of our shale and sandstone, respectively, this represents a generous estimate of the reactive surface area for these rocks.
(d) Erosion and runoff corrections
Many studies have provided evidence that chemical weathering is related to physical erosion [49–54]. Montgomery & Brandon  suggest that erosion is related to relief. Their work in the Olympic Mountains revealed the following relationship: 2.4where E is erosion, E0 is background erosion due to chemical weathering (0.01 mm yr−1), K = 0.00025 mm yr−1 is a rate constant, Rz is the local mean relief (m) and Rc is a limiting value for relief, which we have set to 2000 m rather than 1500 m  to avoid numerical problems with our steepest grid points.
To incorporate the effect of relief on physical erosion, we use the standard deviation of the orography (m) as an approximation of the local mean relief for each individual grid point. Erosion values are then normalized to a mean erosion Emean, calculated by applying equation (2.4) to the mean value of the standard deviation of orography for the world.
Following Millot et al. , we can calculate an approximate erosion correction to our weathering fluxes: 2.5where i refers to the individual grid point.
A further correction is made to account for weathering responses to runoff, which may be linked to effective regolith depth. Berner  combined an empirical relationship between concentration and runoff: 2.6with an approximation for the bicarbonate weathering flux at steady state: 2.7where A is the land area being weathered and V is the volume of water in the regolith, to arrive at the following normalized correction: 2.8where runoff0 and runofft are the global runoffs at 0 Ma and time t, respectively.
Runoff is a major control on basalt weathering . For a given temperature, Dessert et al.  noted that bicarbonate concentrations for a given temperature appear to be constant regardless of runoff, whereas bicarbonate fluxes are proportional to runoff. Although Oliva et al.  report that silica and cation fluxes from granite weathering may rise with increasing runoff where the mean annual temperature is greater than 13°C, concentrations decrease with increasing runoff. No simple relationship between cation fluxes and runoff is apparent for granite , unless physical erosion rates are high . However, Amiotte-Suchet et al.  found that weathering rates of pure silicate sandstones and granites were about five times less responsive to runoff than pure silicate shales and basalts. As we are following the same approach assuming that these rocks comprise exclusively silicate minerals, we apply our runoff correction only to basalt and shale. We adopt the GEOCARB approach to runoff, such that the runoff for each individual grid point is normalized to the maximum runoff for the world. Our correction, therefore, either reduces weathering rates or leaves them unchanged.
By these means, we reduce the weathering contribution from low-relief regions likely to have developed thick soils that protect bedrock from further weathering , and we account for flushing or effective regolith depth variations on basalt and shale. With these corrections, the final weathering rate for each basalt or shale grid point is 2.9where subscripts have the same meanings as for equation (2.3). In summary, the rate constants kij, reaction orders nij and apparent activation energies Eapp,i for weathering agents aj = H+, H2O and OH− operating on mineral i are from the compilation of Palandri & Kharaka , monthly average temperature T is generated for each grid point by HadCM3L, Ce ultimately depends on the orography data, which are boundary conditions for HadCM3L, and Cr depends on the runoff calculated by SDGVM.
3. Model validation
We used modern observational monthly temperature, humidity and precipitation datasets  at a resolution of 1° × 1° to drive SDGVM, producing the NPP map shown in figure 2a. Using the Water Resources eAtlas , we assigned grid points to world watersheds, excluding any watersheds with less than 30 per cent natural vegetation such as forest or grasslands.
Watersheds were considered to be EM if they had at least 50 per cent forest cover comprising tree families such as Pinaceae, Fagaceae and Dipterocarpaceae, which are known to be EM . Outside the watersheds, grid points were treated as EM if the trees were deciduous (such as oaks and larches) or evergreen needle-leaved (such as pines), as shown in table 1.
A lithology map at the same resolution was kindly supplied by Philippe Amiotte-Suchet . We treated all volcanic rocks as basalt; otherwise, we used the rock types of Amiotte-Suchet et al. . In a few cases, we altered the geology to match literature data. For example, the Hong watershed was entirely shale on the lithological map, even though there are igneous and metamorphic rocks in the Thao tributary contributing to the flux of cations in the river . In both the Hong and Amur watersheds, we made adjustments to the lithology so that it would be similar to that of Moon et al. [64,65].
We represent most woody plants as EM except for broad-leaved evergreen trees that most closely match the tropical rainforest distribution map of Read . He considered this biome to be AM with some EM. Our approach will misrepresent EM broad-leaved evergreen dipterocarps and oaks as AM, and could wrongly treat important Late Cretaceous AM needle-leaved trees such as Sequoia  as EM.
(a) Observed and predicted watershed Ca + Mg weathering fluxes
We compared our simulated weathering fluxes of Ca + Mg at the catchment scale with measurements reported in a global compendium  in which 29 watersheds fit the criteria for at least 30 per cent natural vegetation outlined above (figure 3). These watersheds are in Africa, Asia, North and South America and Australia, and contain a range of vegetation types as the dominant land cover, including boreal, temperate and tropical forests. Overall, following the application of erosion and runoff corrections, our predicted fluxes show close correspondence with the measurements reported by Gaillardet et al.  for watersheds producing high Ca + Mg fluxes (R2 = 0.726 for untransformed data). A linear fit comparing modelled and observed log-transformed fluxes has a slope of 0.913; the slope would be 1.0 for a perfect match (figure 3e). The spread of the predictions at the lower end of the scale, however, suggests that the rocks undergoing weathering in these watersheds are not captured by the 1° × 1° lithology map of Amiotte-Suchet et al.  or that our erosion and runoff corrections do not capture smaller scale variations in Ca + Mg flux from silicate weathering. These results are obtained with an input climate averaged over 11 years (1992–2002). We also tested a HadCM3L climate for the modern day that gives a slope of 0.732 for log-transformed data. However, R2 for the untransformed data is only 0.5 for this model run, because the observed Fly River (Papua New Guinea) runoff is over six times larger than that produced by HadCM3L.
These model-data comparisons serve to verify our modelling approach at the catchment scale, despite a number of simplifying assumptions, such as uniform soil type and depth, lithology and mineralogy. It further indicates that our numerical representation of biological and abiotic weathering processes are probably realistic and therefore have some utility when extrapolated into the past. This view is supported by the similarity between our Ca + Mg flux modelled for the Orinoco watershed in South America (0.056 mol Ca + Mg m−2 yr−1) with that obtained from a more detailed multi-layered soil chemistry model (0.045 mol Ca + Mg m−2 yr−1) for the same catchment that includes a treatment of the precipitation of secondary clay minerals, the inhibiting action of aluminium on weathering rates and trace quantities of calcite in silicate rocks .
Because our model simulates seasonal cycles in biology and climate, calculated monthly Ca and Mg weathering fluxes can be compared with observed seasonal fluctuations. In the Hong watershed (figure 4) that lies within China and Vietnam, our simulated seasonal amplitude is similar to observations reported by Moon et al.  showing warmer summers are associated with doubling the weathering fluxes compared with those in the wintertime. It is important to appreciate that modelled catchment scale estimates of cationic weathering fluxes in figure 4 are the mean of multiple grid points with tropical to subtropical climates. Therefore, we explored how our modelled variability compares with observations made under the contrasting conditions of temperate to boreal forest in the Amur watershed in Mongolia . Overall, the comparison shows that the spatial variability of our cationic fluxes is comparable with observations for this watershed.
(b) Global continental weathering
We next extend our simulations of terrestrial NPP and continental Ca + Mg fluxes to the global scale (figure 2b). Comparison of the two global maps indicates a general correspondence between areas of high productivity and high weathering rates, but these rates are also linked to the underlying lithology; regions with the highest weathering fluxes of less than 0.4 mol m−2 yr−1 typically coincide with the location of basalt.
At the global scale, little is known about the spatial extent of highly active weathering regions, making validation of our global weathering map difficult. Hartmann et al.  published a global map of CO2 consumption on a 1° × 1° grid resolution based on applying a spatially explicit function for CO2 consumption by weathering of silicate and carbonate rocks. Our weathering map is not directly comparable with that of Hartmann et al. , which accounts for the weathering of all rock types, including carbonates. Our measure of silicate weathering is the flux of Ca + Mg, which we do not translate into CO2 consumption. Nevertheless, both maps show regions of agreement. In particular, the locations of high CO2 consumption by chemical weathering in the tropical forests of South America, central Africa and throughout South East Asia correspond to our high rates of weathering (figure 2b). Our erosion correction indicates low fluxes for large parts of the northern hemisphere, whereas our runoff correction sharply reduces the flux from the Siberian, Ethiopian and Deccan Traps. We also come to a similar conclusion that weathering hotspots produce the bulk of the global weathering flux. Our results indicate that 10 per cent of the silicate land area of 106 Tm2 in our map produces 66 per cent of the Ca + Mg flux from silicate weathering. This is higher that the figures given by Hartmann et al. , where 9 per cent of the global exorheic area of 121 Tm2 produces 50 per cent of the CO2 consumption. In our case, hotspots occur where relief, runoff and NPP are high, especially on basalt, and it is therefore instructive to investigate the effect of lithology. For an all-granite world of total area 127 Tm2, our model produces 41 per cent of global Ca + Mg from 12.7 Tm2, or 60 per cent of global Ca + Mg from 20 per cent of the land area. Such a world is less influenced by hotspots than expected by Hartmann et al. , and because granites produce more Ca and Mg than sandstones and shales, we get a higher global Ca + Mg flux of 5.20 mol m−2 yr−1. This is closer to some global literature estimates [3,68] but far in excess of estimates for igneous rocks alone (table 2). Furthermore, Dessert et al. [56,57] expected the contribution of basalts to global CO2 consumption to be non-negligible because they weather 5–10 times faster than granite and gneiss. For the 47 grid points on the Deccan Traps, the ratio of the weathering rates for the standard and all-granite runs had a mean of 2.9 and a maximum of 9, with six grid points within the range 5–10 estimated by Dessert et al. [56,57]. Worldwide, our basalt–granite weathering ratio has a mean of 4.9, in line with that reported by Dessert et al. [56,57].
A further approach to validation is possible by comparing our global continental fluxes with earlier published estimates based on elemental ratios and monolithological watershed data to estimate the silicate contribution of solutes to river water (table 2). Global totals are similar whether the model is driven by observations of modern climatology or with the HadCM3L simulated modern climate (table 2). Simulated global weathering of granite and basalt (i.e. igneous) silicate rocks are comparable with estimates published by Hartmann et al.  and Meybeck , after corrections owing to erosion and runoff, as applicable. Strontium isotope data have been interpreted to suggest a present-day basalt flux that is ca 44 per cent  or 35 per cent  of the total, in agreement with Dessert et al.  who estimated the contribution of basalt weathering to global silicate CO2 consumption to be ca 35 per cent. Our results are in agreement with these studies (table 2), with basalt weathering being 45 per cent of the total (table 2).
Our modern flux of 3.77 Tmol Ca + Mg yr−1 implies a CO2 consumption flux owing to the weathering of Ca and Mg silicates of 7.54 Tmol C yr−1 (table 2), which is comparable with the estimate of 7.1 Tmol C yr−1 given by Gaillardet et al. . The CO2 consumption flux from basalt reported by Hartmann et al.  (3.55 Tmol C yr−1) provides an upper limit to the flux of Ca and Mg from basalt of 1.78 Tmol Ca + Mg yr−1, if one assumes that all of the bicarbonate is balanced by Ca and Mg, which again supports our estimate of 1.68 Tmol Ca + Mg yr−1 (table 2).
There are few global estimates of silicate Ca + Mg fluxes from sedimentary rocks (shales and sandstones). Berner et al.  calculated global weathering rates using riverine data of Meybeck  and the assumption that 15 per cent of calcium in all sediments globally is found in silicate minerals . However, this approach gives numbers that are about an order of magnitude larger than our estimated weathering of sedimentary rocks (following erosion and runoff corrections), in spite of our using similar mineralogies for sandstones and shales (table 2). We also employed a very high mineral surface area for these rocks, but the minerals (albite, K-feldspar, chlorite and smectite) have relatively slow weathering rate laws . Without any corrections owing to erosion or runoff, we obtain for these rocks a flux similar to that of Berner et al. . However, we can see no rationale for excluding the effects of erosion and runoff because this would imply that the weathering front on shales and sandstones is never obstructed by thick soils or other occlusions, and that these rocks do not respond to the exposure of fresh surfaces following physical erosion.
4. Present-day weathering by plants and mycorrhizal fungi
As expected from a weathering model driven by the biological proton cycle , our simulations reveal strong relationships between vegetation NPP and the flux of Ca and Mg owing to the weathering of silicate rocks . These relationships vary with lithology and plant and mycorrhizal functional type (figure 5). Because the runoff and erosion corrections (equations (2.5) and (2.8)) account for changes in soil properties such as reacting mineral surface area rather than characteristics of the plants or their symbionts, we analyse the potential (uncorrected) flux of Ca + Mg from our weathering model.
As expected from geochemical theory  and field observations , basalt weathers faster than granite. In the model, this occurs for two reasons. First, because the reactive mineral surface area of basalt is two orders of magnitude higher than granite, and second, because the rate laws for dissolution of minerals of basalt are faster than those of granite (figure 5). However, the higher buffering capacity of basalt compared with granite limits the effect of acidification on its weathering rates. This means that on basalt, acidification by EM fungi is primarily important at low NPP. On granite, the ability of EM fungi to focus protons on minerals makes them more effective weathering agents than AM fungi as NPP increases.
The gradient in NPP at the global scale (figure 2a) is driven primarily by climate linked to vegetation type. Evergreen broad-leaved trees (i.e. tropical rainforests) have the highest NPP of all the simulated PFTs, and they drive the highest weathering fluxes on both rock types (figure 5a,c). Deciduous broad-leaved trees (i.e. temperate forests) are the next most productive, and needle-leaved trees of boreal forests are generally less productive again. In terms of weathering, evergreen broad-leaved trees achieved high rates, even though they are usually AM, because of high productivity driving the biological proton cycle and high acidity of the soil solution.
Deciduous broad-leaved trees often produce lower Ca + Mg fluxes than needle-leaved trees of the same or even lower NPP that occupy cooler climates on basalt. This effect does not occur if all trees are considered to be AM. On granite, all EM trees except for evergreen broad-leaved trees produce Ca + Mg fluxes as a linear function of NPP, and these fluxes are higher than those of AM trees.
The fractional contribution of each plant and mycorrhizal functional type to the global weathered flux of Ca + Mg is given in table 3, which indicates that although they contribute only 46 per cent of global NPP and cover 19 per cent of global land area, evergreen broad-leaved trees contribute 72 per cent of the Ca + Mg flux from silicate weathering. Most of these trees are treated as AM, so a similarly high fraction of the global Ca + Mg flux is due to AM fungi. These results indicate that this PFT is the single most important agent of biological weathering, with mycorrhizal functional type being of secondary importance.
One of the 12 testable hypotheses on the geobiology of weathering  is that the location and extent of biological weathering is controlled by carbon fixation [15,43]. Our model is designed so that biological weathering depends on NPP, but at the global scale, weathering is enhanced at different rates depending on rock type and on plant and mycorrhizal functional types. Comparison of the maps in figure 2 indicates that global NPP maps cannot be used to identify the locations producing the highest silicate weathering fluxes owing to the strong influence of lithology, erosion and runoff. The enhancement of silicate weathering by EM fungi compared with AM fungi is attenuated for the most productive trees on the rocks rich in Ca and Mg. Therefore, figure 5 shows that the hypothesized positive relationship between NPP and biological weathering  is strongly mediated by plant and mycorrhizal functional type acting in concert with abiotic factors.
5. Weathering in the past
Because of the lack of lithological maps for the time intervals corresponding to the suite of HadCM3L palaeoclimate simulations required for continental weathering calculations, we developed our own, starting with entirely granite worlds and adding other rock types by hand. Granite was replaced with flood basalts of appropriate age and size based on published data concerning large igneous provinces [74–78] and websites hosted by the universities of Leicester (http://www.le.ac.uk/gl/ads/SiberianTraps/Index.html) and Bristol (http://palaeo.gly.bris.ac.uk/Palaeofiles/Permian/SiberianTraps.html). Emplacement of the Deccan Traps occurred near the K/T boundary  just after the date of our Maastrichtian run (67 Ma). We have included them in our weathering calculations, and they comprise an important weathering hotspot at this time, and also in the Eocene.
Using the online maps produced by the United States Geological Survey and the Atlas of Mesozoic and Cenozoic Coastlines , we replaced other granitic regions with limestones, shales and sandstones of appropriate age. We attempted to locate sedimentary basins using the available boundary conditions (orography and its standard deviation), seeking simple rules for the probable locations of terrestrially derived sandstones and shales based on their modern distributions . Overall, the lithology of 40–50% of the continental grid points could be predicted by assuming shale for orography under 200 m and standard deviation under 600 m, and sandstone for orography between 200 and 1300 m and standard deviation under 100 m. This approach is distinct from previous attempts to model global continental weathering over geological time, which assumed all grid points contained all rock types in their modern global proportions . However, we aimed to capture the effect of basalt weathering with respect to changes in basalt production and erosion over geological time, as well as palaeogeography and palaeoclimate.
In spite of adding sediments at the expense of granite to generate palaeolithologies, overall land areas underlain by granite were higher than those of Bluth & Kump . However, the latter are known to be underestimated owing both to lack of data about older lithofacies and inability to account for erosion of rock units . Gibbs & Kump  revised the Bluth & Kump modern lithological map, but over one-quarter of their land area is assigned to ‘fold belts’ rather than being included in the lithological categories of Amiotte-Suchet et al. , so their total areas for many of our rock types are lower than ours. We recognize the limitations inherent in this approach to generating lithologial maps for weathering calculations, but the advantages are that it establishes a standardized method for exploring how continental weathering may have changed over geological time. Nevertheless, our estimates of continental weathering should be regarded as exploratory rather than definitive at this stage.
(b) Global patterns of terrestrial vegetation productivity and weathering in the past
Global patterns of terrestrial NPP over geological time vary with latitude, climate, atmospheric CO2 concentration and continental configuration (figure 6). In particular, areas of highest NPP are sustained in the tropics with year-round warmth, high precipitation and high solar energy. Global terrestrial NPP has varied considerably over geological time (figure 6), and functions to account for changes in weathering by plants under different CO2 and climate have been included in models of the long-term carbon cycle [5,8]. Taylor et al.  compared these functions with the estimates that were produced using SDGVM. Total NPP for our time slices compares well with earlier SDGVM estimates of Beerling & Woodward  and with functions developed by Volk  and Bergman et al. , but is considerably higher than that of Berner  when these functions are scaled by modern NPP.
Comparison of global patterns of NPP and weathering in the past serves to reinforce the idea that there is not always a one-to-one correspondence between the two (figures 6 and 7). For example, in the Late Triassic (Norian), although NPP is highest in the evergreen broad-leaved forests around the south-eastern coast of Asia, the highest Ca + Mg fluxes are associated with deciduous broad-leaved trees on the southern Siberian Traps. In the Jurassic (Kimmeridgian), the highest Ca + Mg fluxes occur on the Ferrar (Antarctica and Australia) and Karoo (southern Africa and Antarctica) basaltic regions under herbaceous vegetation, and on the Central Atlantic Magmatic Province (CAMP) basalts under evergreen broad-leaved forests. These sites are not coincident with global maximum NPP associated with evergreen broad-leaved forests on granites of Southeast Asia and Central America. The low runoff simulated over most of the CAMP basalts and the Siberian Traps means that these regions do not contribute as much Ca and Mg to rivers as might have been expected, given their large area and the equatorial position of CAMP.
Other attempts to model weathering in the Triassic have focused on the Carnian and Rhaetian just before and just after the Norian, respectively. Our global weathering pattern (figure 7g) differs substantially from the CO2 consumption maps of Goddéris et al. , where the highest weathering occurs in the equatorial region. Despite high temperatures and adequate moisture to support deciduous broad-leaved trees at the equator, we have no equatorial weathering hotspots in the Norian because most of the rocks are either granite or sedimentary.
Cretaceous simulations indicate that the CAMP basalts produce the highest weathering fluxes in the Aptian (119 Ma), where they coincide with high runoff and high NPP of broad-leaved evergreen forest. This contrasts with the situation in the Cenomanian, where these rocks are associated with low runoff and sparse forest coverage. The weathering map for the Cenomanian bears some similarity to the CO2 consumption maps of Goddéris et al. . Our map shows a greater proportion of flux from Antarctica, which was forested with evergreen needle-leaved trees in general agreement with fossil evidence . Runoff on the CAMP basalts was low and the remainder of equatorial Africa was not underlain by basalt; so equatorial Africa was not a weathering hotspot despite the high NPP and high predicted CO2 consumption in the simulations of Goddéris et al. . Instead, the highest fluxes of the Cenomanian come from AM evergreen broad-leaved forests of the wet Emeishan (China) basalts, followed by EM deciduous broad-leaved and evergreen needle-leaved trees on the somewhat drier and cooler Karoo and Siberian Traps basalts, respectively. EM trees dominate 70 per cent of the basalt grid points in the Cenomanian, in agreement with Taylor et al. , but this is not the case in the latest Maastrichtian when the CAMP and Deccan Trap basalts are covered with high-NPP AM evergreen broad-leaved trees.
(c) Biological weathering on the continents over geological time
For each time step, we undertook three simulations to enable us to investigate and isolate the changing role of vegetation, mycorrhiza and climate on continental weathering rates over geological time: (i) a standard run with mycorrhizal functional type allocated to trees according to table 1 (all AM in the Aptian, Kimmeridgian and Norian); (ii) a non-mycorrhizal run, with only fine roots directly impacting soil chemistry; and (iii) a run with no impact of roots or fungal hyphae on soil chemistry. The latter ‘no vegetation’ run still includes the soil water content and runoff from SDGVM. The net direct effect of land plants on chemical weathering is calculated as the difference between the total global Ca + Mg fluxes for the mycorrhizal and ‘no vegetation’ cases. Results, along with the atmospheric CO2 used to force the models, are shown in table 4.
Figure 8a indicates the potential for up to fourfold variations in continental weathering over the past 220 Myr relative to the present day (table 4). These absolute global fluxes obtained from our process-based modelling are of the same order of magnitude as those obtained with global geochemical models using independent methodologies [8,9,29], but are often higher, as shown in figure 8a. According to the above sensitivity analyses, climate dominates global weathering reflecting the strong temperature control on rate laws of mineral dissolution (figure 8b,c). The exceptions to this pattern are the Late Triassic (Norian, 216 Ma) and Early Jurassic (Kimmeridgian, 153 Ma) simulations, which show unusually low weathering regimes, given the warm average land surface temperatures. In these two cases, the continental configuration of land masses created a supercontinent with extensive deserts and a more arid environment (figure 8c), causing low soil water content and runoff that hinders chemical weathering.
According to our modelling, vegetation enhances weathering over and above that caused by climate in all simulations, in agreement with numerous field studies [15,83] (figure 8a,b). For the contemporary climate and CO2 (i.e. 0 Ma), mycorrhizal vegetation increases continental weathering by a factor of 1.3 at the global scale, and the same set of simulations afford a further opportunity for independent validation against observations. On Icelandic basalt at 0 Ma, this ratio increases to 2.8 for sites dominated by deciduous broad-leaved trees, and 1.9 for those with evergreen needle-leaved trees. These numbers are slightly low when compared with estimates of three to four derived for adjacent-vegetated and non-vegetated catchments in Iceland . Moulton et al.  did not find a significant difference between deciduous broad-leaved trees (Betula pubescens) and evergreen needle-leaved trees (Picea spp., Abies spp.), although they note that the B. pubescens trees were stunted.
In the Norian (216 Ma), Kimmeridgian (153 Ma) and Aptian (119 Ma), mycorrhizal fungi exert a minor additional effect on biological weathering because they are exclusively modelled to be the ancestral AM type that predominated prior to the later evolution of EM fungi. AM fungi are thought to be less effective than EM fungi because they enhance rather than replace the nutrient uptake functions of fine roots, and because they do not tend to exude organic acids .
In all simulations, after and including the Cenomanian, the boreal and temperate forests are assigned EM fungal partnerships following the function of Taylor et al.  that characterized the sigmoidal rise in EM vegetation at weathering hotspots during the Cretaceous . This function has 0.1 per cent EM vegetation at weathering hotspots at the estimated time of EM origination (180 Ma) according to molecular data [15,85], rising to 71 per cent at 95 Ma. As a result of the ability of these fungi to target nutrient-bearing minerals, EM and AM fungi together contribute as much as or more than the amount of weathering driven by vegetation alone (figure 8b). In the Cenomanian, mycorrhizal fungi exert an especially large effect (figure 8a,b), and EM vegetation is found on 70 per cent of basalt grid points, including the Karoo and Siberian Traps as described above, of which more than 40 per cent is provided owing to fungal activities (figure 8b and table 4). These global-scale simulations provide new support for the idea previously proposed by Taylor et al. : the rise of EM trees may have contributed to the CO2 drawdown previously attributed to the spread of angiosperms since the Mid- to Late Cretaceous. From these results, we can create a new dimensionless plant weathering efficiency function (fE) to better characterize how terrestrial ecosystem weathering changed over geological time (figure 8d). The resulting pattern shows that terrestrial ecosystem weathering could have been up to four to five times more efficient in the past relative to today, especially during times of global warmth, high CO2 and extensive land cover by EM vegetation. The pattern of ecosystem weathering efficiency over the past 215 Myr differs substantially from the fE function included in geochemical models [7,8], which is between 0.875 and 1.0 for the past 250 Myr. Because fE is normalized to the present-day flux, our values exceed one in the Cretaceous. One possible explanation for such unexpectedly strong plant effects is that the CO2 fertilization function of the GEOCARB models does not adequately account for the significantly higher global NPP of the Cretaceous. Cenomanian global NPP (138 Pg C yr−1) is more than twice as high as it is today (62 Pg C yr−1).
6. Concluding remarks
Our global modelling strategy is well-validated against modern observations at the catchment and global scales and, therefore, appears to be robust, at least for the present-day situation (figures 3,4 and table 2). For contemporary ecosystems, our simulations indicate that the mycorrhizal status of trees has an effect on spatially resolved Ca + Mg fluxes from silicate weathering. The more ‘aggressive’ EM fungi do not necessarily produce higher fluxes if the trees have high NPP and are weathering basalt (figure 5). For NPP under 10 Mg C ha−1 yr−1, EM fungi enhance weathering compared to AM fungi (figure 5). Fine roots and hyphae enhance weathering fluxes of Ca and Mg by a factor of 1.2–2.1, and the highest fluxes occur where NPP is the highest for a given rock type. Mycorrhizal functional type is however less important to global weathering fluxes than PFT in the modern world.
When extrapolated back into the past, the same modelling approach indicates the potential for major enhancement of weathering by terrestrial ecosystems, with the coevolution of ectomycorrhizas and vegetation being an especially important factor (figure 8 and table 4). The enhancement is greater than previously suggested by empirical functions  with implications for the role of the evolving biota in regulating the global atmospheric CO2 concentration. One of the major differences between this study and previous studies of the effect of plants on the long-term geochemical carbon cycle [7–10,20] is that we present fluxes (tables 2–4; figures 2b, 3–5, 7 and 8) that have been calculated without calibration to any assumed or observed weathering rate or degassing flux. Our model could be calibrated by changing the reactive mineral surface area for each rock type until our modern flux matched that of Gaillardet et al. , for example. Such a calibration is probably necessary if we were to use our model to generate revised estimates of CO2 drawdown to produce a new CO2 curve, because our uncalibrated fluxes could cause unrealistic CO2 drawdown.
In conclusion, this study demonstrates ‘proof of concept’ by reporting novel approaches for the more realistic process-based treatment of plant–fungal interactions with the geosphere at the global scale. Accurately capturing such interactions is critical for better understanding the coevolution of terrestrial ecosystems, CO2 and climate over geological time, and constitutes a first step towards developing ‘next-generation’ geochemical models incorporating the effects of mycorrhizal evolution.
We thank Philippe Amiotte-Suchet for providing his lithological map, and Mark Lomas for documentation related to SDGVM. We acknowledge funding of this research through NERC award NE/E015190/1 and NERC/World Universities Network consortium award NE/C521001/1. Lyla Taylor was supported by a Hossein-Farmy studentship at the University of Sheffield linked to these grants. David Beerling gratefully acknowledges support from a Royal Society–Wolfson Research Merit Award. We are grateful to Robert Berner and Yves Goddéris for helpful reviews of our manuscript.
One contribution of 12 to a Theme Issue ‘Atmospheric CO2 and the evolution of photosynthetic eukaryotes: from enzymes to ecosystems’.
- This journal is © 2012 The Royal Society