The Rock Geochemical Model ( RokGeM ) v 0 . 9

A new model of terrestrial rock weathering – the Rock Geochemical Model (RokGeM) – was developed for incorporation into the GENIE Earth System modelling framework. In this paper we describe the model. We consider a range of previously devised parameterizations, ranging from simple dependencies on global mean temperature followingBerner et al. (1983), to spatially explicit dependencies on run-off and temperature (GKWM, Bluth and Kump, 1994; GEM-CO2,Amiotte-Suchet et al. , 2003) – fields provided by the energy-moisture balance atmosphere model component in GENIE. Using long-term carbon cycle perturbation experiments, we test the effects of a wide range of model parameters, including whether or not the atmosphere was “shortcircuited” in the carbon cycle; the sensitivity and feedback strength of temperature and run-off on carbonate and silicate weathering; different river-routing schemes; 0-D (global average) vs. 2-D (spatially explicit) weathering schemes; and the lithology dependence of weathering. Included are details of how to run the model and visualize the results.


Motivation
A new model of weathering is developed for incorporation into the GENIE Earth System Model in order to allow for the simulation of spatially explicit, full global carbon cycle perturbations to complete equilibrium on a million year timescale.Our aim was to produce an intermediate complexity model of weathering, which advances on 0-D global models such as GEOCARB (Berner, 1990(Berner, , 1991)), but can still be integrated for a million years (to equilibrium) and is thus less complex than process-based models such as B-WITCH (Roelandt et al., 2010;Beaulieu et al., 2012).
The carbon cycle involves short-term (years to centuries) terrestrial processes, medium-term (decades to millennia) oceanic processes, and long-term (10 3 -10 5 yr) geologic processes.The latter involve the weathering of carbonate and silicate rocks.Carbon dioxide dissolves in rainwater to form a weak carbonic acid that breaks down carbonate and silicate rocks, producing calcium, magnesium and bicarbonate ions.Through riverine transport, this bicarbonate reaches the ocean.There the calcium and magnesium ions recombine with bicarbonate to form new carbonate sediments (Ridgwell and Edwards, 2007).In the case of carbonate weathering there is an overall null cycle for CO 2 , whereas silicate weathering transfers CO 2 to the Earth's crust.
Examples of simulations possible with the new (GENIE-)RokGeM model include the eventual return to pre-industrial concentrations of atmospheric pCO 2 (over the course of 1 Myr) following anthropogenic emissions, and Phanerozoic events such as the Paleocene-Eocene Thermal Maximum.Results from fossil fuel perturbations are presented here.

Early modelling: the carbon cycle over Earth history
Early work using box models (not spatially resolved) now known as WHAK (Walker et al., 1981), BLAG (Berner et al., 1983) and GEOCARB (Berner, 1990(Berner, , 1991) ) demonstrated the essential role silicate weathering plays in the regulation of climate over geological timescales.In WHAK, the global rate of silicate weathering is given a dependence on atmospheric pCO 2 (CO 2 partial pressure) and temperature.The draw-down of atmospheric CO 2 through silicate weathering, and the positive relationship between temperature and pCO 2 , creates the now commonly known silicate weathering negative feedback loop.WHAK employs fairly rudimentary parameters, relying on results from early general circulation models (GCMs), but it is the archetype for models linking global silicate weathering and climate.Following WHAK, it was suggested that the abiological model be amended to include biology, noting that pCO 2 in soil is enhanced to levels much higher (10-40 times) than atmospheric levels due to the presence of terrestrial life (Lovelock and Whitfield, 1982;Lovelock and Watson, 1982).This gives the geological silicate weathering feedback a geophysiological (or Gaian) flavour; it is argued that this is necessary to keep pCO 2 and hence temperature low enough for continued functioning of the biosphere over the eons.It was later calculated, using a model based on WHAK, but including an (estimated) biological amplification factor of 100 to 1000, that an abiotic Earth would be between 30 or 45 • C warmer (Schwartzman and Volk, 1989).The biological amplification factor included the raising of soil pCO 2 as well as the stabilisation of soil; the production of organic acids (a soil matrix has high surface area, which stores water, containing the acids produced in solution which weather the minerals in contact with them); and enhanced physical weathering through the fracturing (roots) and microfracturing (microbes) of rocks and mineral grains.Further work (Schwartzman and Volk, 1991) demonstrated the plausibility of a biologically enhanced silicate weathering feedback being a mechanism for planetary homeostasis over the full range of Earth history.Indeed, applying a similar model eons into the future, the remaining lifespan of the biosphere was calculated to be 0.9-1.5 Gyr (Caldeira and Kasting, 1992).Later modelling (Lenton and von Bloh, 2001), including a direct (hyperbolic) dependence of weathering on (biological) productivity, showed that biotic and abiotic states of the Earth system were bistable; with a large enough perturbation (the size of which diminishes with increasing solar luminosity), the system can flip from the biotic to abiotic state.
BLAG put the silicate weathering mediated geological carbon cycle into a late Phanerozoic (last 100 Myr) context.Weathering processes were, as with WHAK, abiological, but precipitation of carbonate in the ocean was included, which is biological in origin.Also included were geochemical representations of the tectonic processes of metamorphic decar-bonation and CO 2 outgassing.The model involves a series of non-linear mass balance equations for determining fluxes into and out of land and ocean mineralogical reservoirs, as well as the atmospheric CO 2 reservoir.BLAG was the basis for the more extensive GEOCARB model, which was constructed to model atmospheric pCO 2 over the complete Phanerozoic eon.GEOCARB simplified the maths of BLAG by aggregating carbonate minerals and combining the ocean and atmosphere into a single reservoir.The soil-biological enhancement of weathering was included as a parameterization.Other factors controlling weathering rates in GEO-CARB are continental land area, mean elevation and runoff.Results of the model show that both geological and biological effects are important in the regulation of CO 2 over geologic time.Successive versions, GEOCARB II (Berner, 1994) and GEOCARB III (Berner and Kothavala, 2001) have refined the model by constraining and detailing its parameterizations with new scientific observations and (for climate variables) GCM output.A larger Phanerozoic geochemical cycling model that includes other elements than carbon, COPSE (Bergman et al., 2004), expands on the weathering feedbacks of GEOCARB with an interactive biota that is coevolutionary with CO 2 .
Moving from the geological, to mere glacial-interglacial timescales, an attempt was made (Munhoven, 2002) to account for the oscillations in atmospheric CO 2 recorded in the Vostok ice cores by applying the Gibbs and Kump Weathering Model (GKWM; Bluth and Kump, 1994) and Global Erosion Model for CO 2 fluxes (GEM-CO2; Amiotte Suchet and Probst, 1995) to GCM output of run-off from the PMIP (Paleo Model Intercomparison Project).It was found that pCO 2 concentrations were drawn down by a small (4.4-13.8ppm), but non-negligible fraction of the glacialinterglacial variation (∼ 75 ppm).However, it is unclear whether there has been significant variation of weathering on glacial-interglacial timescales (Foster and Vance, 2006).GKWM and GEM-CO2 are incorporated into the RokGeM model described in this paper (see Sect. 3.6).

Later modelling: anthropogenic and other Phanerozoic perturbations of the carbon cycle
Various attempts have been made to quantify the lifetime of the anthropogenic CO 2 perturbation.However, the timescale over which the silicate weathering feedback operates has only been quantified once, using a box model, (Sundquist, 1991).The model had a one-box atmosphere, coupled to a mixed layer ocean underlain by an 11-box-deep ocean; each deep ocean box was coupled to a sediment box.Carbon cycle chemistry was performed using interactive variables for atmospheric CO 2 , ocean DIC and alkalinity, and sedimentary CaCO 3 , and the inclusion of simplified equations of carbonate and silicate weathering.Expressed as an e-folding timescale, the silicate weathering timescale was calculated to be in the range 3-4 × 10 5 yr.Other future projections of the lifetime of anthropogenic CO 2 have assumed values of 200 kyr (Archer et al., 1997) and 400 kyr (Archer, 2005).
Taking the parameterizations of weathering temperature dependence from GEOCARB, modified by a plant productivity from COPSE, and applying them to a box model compartmentalised into atmosphere, vegetation, soil and multibox ocean and sediments (Lenton and Britton, 2006), the silicate weathering timescale was determined to be of the order of a million years.The silicate weathering process has yet to have its timescale quantified using a spatially explicit model that includes a 3-D ocean, and climate feedback on temperature and circulation.
Recent work using spatially resolved Earth System Models (ESMs) (Montenegro et al., 2007;Ridgwell and Hargreaves, 2007) has incorporated millennial timescale processes into modelling of the Anthropocene carbon excursion.Detailed, spatially explicit models of the carbonate sediments in the ocean are included, but terrestrial weathering processes are still only dealt with (in the most basic way) as a global average prescribed flux.As models were "only" integrated over timescales in the 10 4 -10 5 yr range, silicate weathering is ignored altogether; the weathering flux is from carbonates, and used to quantify the effect of neutralising fossil fuel carbonic acidity.
An example future projection using one of these models (GENIE-1, the ESM used for the work in this paper, see Sect.1.4 below) has 4000 GtC of emissions sequestered as follows (Ridgwell and Hargreaves, 2007): 66 % are absorbed by the ocean on a timescale of hundreds of years; the dissolution of sea-floor sediments sequesters another 12 % on a timescale of ∼ 1.7 kyr; carbonate weathering on land consumes a further 14 % on a timescale of ∼ 8.3 kyr.Overall, 8 % is left to be removed through silicate weathering.It is worth noting that the proportions of carbon removed from the atmosphere by each of the processes outlined above are different for different sized fossil fuel burns (emissions).For example, with a smaller burn, the proportion of carbon removed by the faster-acting process of ocean invasion is greater, due to the fixed ocean reservoir of alkalinity being relatively larger in comparison to the burn size.
Building on earlier box modelling work (Sect.1.2; this section), and using parameterizations constrained by field data, including the effect of lithology (Gibbs et al., 1999;Amiotte-Suchet et al., 2003), a new weathering model was developed for this paper.This weathering model is incorporated into the GENIE-1 Earth System Model.

The GENIE model and GENIE-1
GENIE -the Grid ENabled Integrated Earth system modelis a modularised framework for Earth system modelling.The constituent components of the Earth system (oceans, atmosphere, land, ocean biogeochemistry, ocean sediments etc.) are encapsulated in separate modules, which can be "plugged and played" together, the idea being that for different uses, different modules can be used, and models of differing levels of complexity can be created.
GENIE belongs to the class of models known as EMICs (Earth system Models of Intermediate Complexity).In contrast to GCMs that use underlying low-level physics equations, EMICs are often highly parameterized.Processes that take place over small spatial and temporal scales are aggregated into high-level parameterizations in order to minimise computational cost.This has the advantage of allowing more processes to be modelled, and also longer integrations and larger ensembles.However, the increase in parameterization has the cost of an increase in the uncertainty of model results.In order to address this issue, during the course of its development, the GENIE model has been tuned using a number of different techniques: an ensemble Kalman filter (Hargreaves et al., 2004); Kriging, response surface models and non-dominated sorting genetic algorithms (Price et al., 2006(Price et al., , 2007)); "pre-calibration" by emulation and plausibility checking (Holden et al., 2009;Edwards et al., 2010); and using Iterative Importance Sampling (Annan and Hargreaves, 2010).
GENIE-1 refers to the first instance of the model created, with a low degree of complexity, and a fast run time.Its basis is a 3-D frictional-geostrophic ocean and an energy moisture balance atmosphere (see below).A thousand years of simulation are capable of being run in one hour on a "normal" desktop PC.For this paper a version of GENIE-1, genie_eb_go_gs_ac_bg_sg_rg, is used.This includes an EMBM (energy moisture balance model) 2-D atmosphere (Weaver et al., 2001) (eb), an 8-layer version of the Goldstein frictional-geostrophic ocean (Edwards and Marsh (2005); and references therein) (go) and sea-ice (gs); AtChem, an atmospheric chemistry module to pass gas fluxes (ac); Bio-GeM (BioGeochemical Model), the ocean biogeochemistry module (Ridgwell et al., 2007) (bg); SedGeM (Sedimentary-Geochemical Model), the ocean sediments model (Ridgwell and Hargreaves, 2007) (sg); and RokGeM (Rock-Geochemical Model) -the model described here (rg).A schematic of the model is shown in Fig. 1.
The "Grid ENabled" part of the acronym GENIE refers to design aspects allowing for the execution of many simultaneous instances of the model across computer clusters, and for the collation, storage, accessing and sharing of results, amongst a geographically spread out community of researchers.Over the course of GENIE's history, a number  of bespoke portals, middlewares, and stand-alone applications/scripts have been developed to enable more efficient use of the modelling framework.Unfortunately, common standards for usage have not been adopted across the entire GENIE community, leading to many different individual efforts at modelling portals, including much work done for this paper (see Sect. 2.3).

Model and experimental set-up
Here the computer model set-up is presented, incorporating its development (Sect.2.1), configuration (Sect.2.2) and running (Sect.2.3).The design and execution of long-term future carbon-cycle perturbation experiments involves a twostage model spin-up (Sect.2.4).Finally, visualisation and some analysis of the results of said experiments is briefly discussed (Sect.2.5).

Model development
The Following software development best practice, the versioning system Subversion was used.Versioning systems allow changes to files to be tracked, and earlier versions to be reverted to at any time, as well as the branching and merging of files, and directory structures.They are indispensable for large software projects, including the development of models containing large code bases, such as GENIE.The Subversion repository for the GENIE project is held on a secure server at the University of Bristol (https: //svn.ggy.bris.ac.uk/websvn/listing.php?repname=genie; username=genie-user, password=g3n1e-user).The main "trunk" of the code base was "branched", creating greg-s_branch for the purposes of developing RokGeM.Regular commits of the code were performed over the course of development to provide both a record, and a means for reversion to earlier versions of code.After significant developments, the development branch was merged back into the trunk, allowing other GENIE users and developers to use the RokGeM model.Also, the trunk was merged back into the development branch at regular intervals in order to avoid incompatibilities causing conflict.In this respect, the suite of tests developed for GENIE were also run before each committal of code to the repository, to minimise the risk of bugs finding their way into the model.The model runs presented in this paper are run using the versions of GENIE contained in greg-s_branch Revs 5439-6324 https://svn.ggy.bris.ac.uk/websvn/listing.php?repname=genie&path=%2Fbranches%2Fgreg-s_ branch%2F&rev=6324.There is minimal change to model output across this span of versions.A tagged releasetags/Greg_Colbourn_PhD_thesis -https://svn.ggy.bris.ac.uk/websvn/listing.php?repname=genie&path= %2Ftags%2FGreg_Colbourn_PhD_thesis%2Fis provided in the svn repository for posterity.
A key input to the weathering model is run-off (see Sect. 3.3,Sect. 3.6).Moisture is transported into continental interiors by the EMBM atmosphere by a combination of diffusion and advection according to prescribed NCEP wind fields.With the values of moisture diffusion and advection given, diffusion is dominant, leading to overly uniform precipitation and hence run-off (see Sect. 3.4 of Colbourn (2011) for details).
Following Ridgwell et al. (2007) there is no seasonality in the model.Climate sensitivity (at equilibrium for a doubling of atmospheric pCO 2 ) in terms of radiative forcing was adjusted from a default of 4 W m −2 (Edwards and Marsh, 2005;Hargreaves et al., 2004;Ridgwell and Hargreaves, 2007;Ridgwell et al., 2007) to the IPCC's estimate of 3.7 W m −2 (Sect.2.3.1 of Solomon et al. (2007); this was done by setting the parameter ea_delf2x to 5.34).

Running the model
The long timescales of the model runs (of the order 10 kyr-1 Myr) presented here made setting off runs and collating results non-trivial.Bash scripts were employed to break runs into manageable pieces, which were submitted as jobs lasting ∼ 5 h each on the ESCluster supercomputer at UEA (see Appendix A1.3 for details).In order to further automate running the model, scripts were also used to automatically spinup the model and set experiments going.A final level of   automation was to generate ensembles and groups of ensembles using scripts written in Mathematica.Results were collated by amending saving routines in the original Fortran code of the relevant GENIE modules (BioGeM, SedGeM and RokGeM) so that time series output was appended to previously existing output files for runs restarted as a result of being divided into manageable pieces.Mathematica scripts were written to collate and visualise data from ensembles; see Appendix A1.6 for details.

Spin-up
The model was spun-up using a two stage process (the technical aspects of which are detailed in cGENIE.HOWTO.pdfSect.8.1, available at http://bit.ly/I4iPJP).Firstly the sediments were "closed" to allow the faster parts of the system (ocean and atmosphere) to equilibrate.The weathering flux is subtracted from ocean cells overlying the sediments to balance the global budget and ensure a closed system.This subtraction involves partitioning the total global weathering flux between each ocean floor cell with a subtraction in proportion to the estimated CaCO 3 preservation and burial rate.In effect, any carbonate hitting the ocean floor is immediately remineralised.For the global average weathering case, terrestrial weathering flux was divided evenly between carbonate and silicate weathering, and set equal to a first approximation of the burial flux of carbonate material in the sediments (F CaSiO 3 ,0 = F CaCO 3 ,0 = 5 Tmol yr −1 ).For the spatially explicit case, weathering fluxes were left to freely equilibrate with the ocean-atmosphere-sediment system.Silicate weathering was balanced by setting outgassing equal to it (F outgas = F CaSiO 3 ,0 ; Eq. 4).The ocean, atmosphere and biogeochemistry were left to equilibrate under the condition of fixing atmospheric pCO 2 at a pre-industrial level (278 ppm).After 25 kyr (ample time for equilibrium to be achieved), the sediments were opened and left for 100 kyr to equilibrate with the rest of the system.For this second part to the spinup, in the global average weathering case, terrestrial weathering flux was set to equal the burial flux of carbonate as diagnosed in the first stage of the spin-up.The system is very nearly, but not quite at equilibrium after this, as shown in Fig. 2, which shows wt % CaCO 3 progressing toward equilibria over the course of 100 kyr second-stage spin-ups.The maximum change in wt % CaCO 3 for any run was less than 0.25 % in the last 10 kyr (after ∼ 75 kyr, it is less than 1 %, so future second-part spin-ups will be kept to this length).Under control conditions (no emissions), the model remained in a stable pre-industrial state after turning off the restorative forcing used in the spin-up.
Where set, the weathering feedbacks had very little effect on the outcome of the spin-up, on account of baseline parameters being set to their global averages (as diagnosed in the first part) in the second part.
It was thought that a shortcut in running experiments might be introduced by way of having a pre-prepared baseline spinup available for use with different ensembles of experiments.However, it was found that, despite outputs being similar in terms of state after spin-up for most model runs (the large range in Fig. 2 is mostly down to one particular ensemble), differing settings for parameters between spin-ups and experimental runs makes for unstable experimental runs, which do Geosci.Model Dev., 6, 1543-1573, 2013 www.geosci-model-dev.net/6/1543/2013/not tend toward equilibrium.Even different runs in the same ensemble have slightly different outcomes when being run from a single spin-up rather than separate spin-ups (a test showed a difference of 0.07 % after 100 kyr when looking at atmospheric pCO 2 ).

Visualising results
With the production of vast amounts of data from model runs numbering over a hundred, displaying results in a systematic, coherent and concise manner is a challenge.In the following section (Sect.3), an overview of results is shown in the form of a plot of pCO 2 for each section describing a model variable.Tables of key numbers from graphs are also produced, presented in Appendix B (Supplement).
The nature of the experiment, involving a restoration of the Earth System to a preindustrial climate state from pulses of carbon emissions, calls for pertinent data such as percentages of excess CO 2 as well as absolute values of CO 2 at specific times, to be tabulated.Also tabulated are the inverse forms -years where specific target values and fractions remaining are reached.These data are also tabulated for surface global warming, and surface ocean acidification, parameters of wide societal (as well as ecological) significance.A full set of results are in Appendix B (Supplement).
On a spatial level, gigabytes of NetCDF data-sets were produced as part of the model output.These were converted into time-varying video visualisations including latitudinal trend-lines and comparison amongst ensemble members.A complete set of these animations is available on request.For the electronic version of this paper, links to selected animations are provided at the relevant places.

Time series
The nature of the experiments conducted for this paper lead to results that are non-linear over time.Carbon cycle perturbations from large emissions pulses are non-linear because there are processes involved that act at different rates.Ocean invasion and mixing are relatively quick (10 0 -10 3 yr timescales), whereas ocean sediment dissolution and terrestrial weathering are relatively slow (10 3 -10 5 yr timescales).In order to present line graphs that show in detail the action that is happening, a number of different approaches were tried.These are presented in Fig. 3. Plots with a linear timescale (Fig. 3a) are dominated by the "long tail" (Archer et al., 2009) of the perturbation.Using multiple timescales in a broken-up plot (Fig. 3c) is a good way to visualise processes happening on different timescales, without resorting to a logarithmic timescale (Fig. 3b), which is less intuitive to grasp (in everyday thought we are accustomed to thinking of time in a linear way).For the dependent variable, a logarithmic scale (Fig. 3d, e, f), or even multiple linear scales (Fig. 3g), can also be used to show more detail.Alternately, separate plots can be shown for each emissions scenario (Fig. 3h), allowing detail to be seen for the much smaller 1000 GtC scenario.
It was decided on balance that using a different linear scale for each emissions scenario for the dependent variable and multiple linear scales for the time axes (Fig. 3h) allowed plots to provide the most information in an intuitive manner, whilst retaining a somewhat uncluttered look.All further time series results are presented in this manner.

Carbonate and silicate weathering
Long-term sequestering of atmospheric CO 2 occurs through the weathering of terrestrial (calcium-or magnesium-) carbonate and silicate rocks.Firstly, in the case of carbonate rock weathering, carbon dioxide dissolves in rainwater to form a weak carbonic acid that breaks down the rocks, producing calcium and bicarbonate ions (Eq.1; here given for calcium, but could equally be for magnesium).Through riverine transport, this bicarbonate reaches the ocean (Archer et al., 1997;Ridgwell and Edwards, 2007).
The partitioning of carbon between the atmosphere, ocean and sediments still leaves a small fraction of initial emissions remaining in the atmosphere after equilibrium is reached on a timescale of around 100 kyr (Goodwin and Ridgwell, 2010).This remaining fraction of anthropogenic carbon is removed from the atmosphere by the process of silicate weathering, whereby two moles of carbon are consumed for every one that is available for subsequent transfer back to the atmosphere from the ocean-sediment system (Eq.2; compare with Eq. 1): Equation ( 2 plagioclase feldspar endmember Anorthite: In the initial zero-dimensional (0-D) instance of the model (Globavg), global average fluxes of calcium ions (Ca 2+ ) from calcium-carbonate (F CaCO 3 ) and calciumsilicate (F CaSiO 3 ) weathering are prescribed.These are modulated by a simple temperature (T), run-off (R) and productivity (P) feedbacks (see Sects.3.2-3.4below) that can be switched on and off.The fluxes are then used to calculate fluxes of DIC (F DIC ) and Alkalinity (F Alk ): Note that there is only one mole of DIC for each mole of Ca 2+ ; this is a short-circuiting of the atmosphere based on the assumption that the atmosphere and surface ocean are well mixed on the timescales considered here.Instead of removing one mole of CO 2 from the atmosphere -and by implication the ocean -and adding two moles of bicarbonate to the ocean (as in Eq. 1), nothing is taken from the atmosphere and one mole of bicarbonate is added to the ocean.The same short-circuiting leaves no net DIC flux for silicate weathering (the DIC terms cancel in Eq. 2).F outgas is a prescribed input of carbon from volcanic outgassing used to counter the prescribed silicate weathering flux (which removes carbon into the geologic reservoir) so as to keep the system in long-term equilibrium.
As shown in Fig. 4, short-circuiting the atmosphere has little effect on long-term global atmospheric pCO 2 levels.Even on shortest timescales, there is minimal difference between the runs with short-circuiting on, and those with it off.However, it is interesting to note that there is a significant difference in alkalinity flux to the ocean (∼ 10 %) between the www.geosci-model-dev.net/6/1543/2013/Geosci.Model Dev., 6, 1543-1573, 2013 two (see Supplement, Fig. 23e), which leads to a significant difference in sediment preservation (Fig. 23d).

Temperature dependence of weathering
Temperature affects weathering through controlling the rate of the chemical reactions involved.Following Berner's geologic carbon cycle models (such as BLAG, Berner et al., 1983 andGEOCARB, Berner, 1994), the dependence of carbonate weathering on temperature is given by where T is temperature; 0 stands for the initial value, and the constant k Ca = 0.049.This equation is derived by correlating concentrations of bicarbonate in groundwater with the temperature of the groundwater (Harmon et al., 1975).Laboratory studies of the temperature dependence of the dissolution of Ca and Mg silicates form the basis of the silicate weathering-temperature feedback used (Brady, 1991): which is an Arrhenius rate law equation; R is the molar gas constant (not to be confused with run-off as it is in all other uses apart from this section) and E a the activation energy for dissolution.Assuming that T ≈ T 0 , this equation simplifies to where for E a given in kJ mol −1 and T 0 in • K.The error induced by this simplification (which lessens the computational resources used by the model -the factor k T can be calculated at initialisation), is less than 10 % for −26 • C < (T − T 0 ) < 31 • C (see Fig. 5).A best guess of E a = 63 kJ mol −1 (Brady, 1991) is used, along with T 0 = 288 K (global average preindustrial temperature), to give k T = 0.09 for the model runs presented in this paper.Another estimate, based on least squares fitting of field based weathering data to multiple controlling variables (temperature, rainfall and erosion), is E a = 74 ± 29 kJ mol −1 (West et al., 2005).The error gives this estimate of E a a large range, spanning more than a factor of two, giving very different weathering flux estimates when the exponential function is taken into account, as shown in Fig. 5.The effect of the temperature feedback is shown in Fig. 6.The silicate weathering-temperature feedback (T_Si) has a greater effect on the draw-down of atmospheric pCO 2 than the carbonate weathering-temperature feedback (T_Ca).The removal of CO 2 from the atmosphere by way of carbonate weathering is completed after ca.50 kyr; this is when the carbonate ocean sediments are restored, and the burial rate of carbon matches the riverine input of carbon to the ocean from the weathering of carbonate rocks.The case with T_Ca switched on leads to a very slight reduction in the equilibrium pCO 2 (312 ppm vs. 314 ppm) when compared to the case with no feedback.
With both weathering-temperature feedbacks switched on, a sensitivity test was conducted over the range of different estimates for activation energy (E a ).As shown in Fig. 7, there is a significant difference between the ensemble members in atmospheric pCO2 drawdown by year 100 k; a range of 50 ppm (13 %) for the 5000 GtC scenario, compared with 81 ppm (18 %) as the difference between the runs with and without weathering-temperature feedbacks (Fig. 6).

Run-off dependence of weathering (0-D scheme)
Weathering is dependent on the residence time of water in contact with rocks.Static water in contact with rocks becomes saturated with solutes; further weathering occurs when water is removed and replenished.Run-off can be used as a proxy for this water cycling that is key to rock weathering.Run-off dependence is factored in following Berner (1994): where R is run-off.The difference between the formulations for carbonate and silicate weathering stems from the assumption that carbonates saturate ground water in contact with them, so that concentrations of weathered bicarbonate are not assumed to be dependent upon run-off (Berner et al., 1983), whereas concentrations of bicarbonate from the weathering of silicate rocks are assumed to be diluted by runoff; β is taken to be 0.65 (Berner, 1994) as a default for the model runs presented in this paper.More recent work gives β = 0.80±0.32(West et al., 2005).The effect of changing β is modest.As illustrated in Fig. 8, there is an ∼ 85 % difference in feedback strength between the end members at 25 • C below the baseline temperature (T 0 ), and only a ∼ 35 % difference at 25 • C above T 0 .Run-off is taken explicitly from the coupled EMBM; or as an alternative option, can be parameterized as a function of temperature following the approach of Berner:   7) (solid lines) to Eq. (8) (dashed lines).Values for E a are taken from (West et al., 2005) (74 ± 29 kJ mol −1 ) and (Brady, 1991) (63 kJ mol −1 ).Note that lower activation energies give bigger fluxes for T < T 0 (bottom panel, different y axis scale).
observations.Figure 9 shows a plot of Eqs. ( 10) and ( 11) with the substitution of Eq. ( 12).Even for extreme values of k run , the feedbacks are much less sensitive than the silicate weathering-temperature feedback is to E a (Fig. 5).
The effect of the run-off feedback is shown in Fig. 10.Only when the emissions pulse is large (5000 GtC) is there a significant effect.Even so, the effect is much smaller than that due to temperature feedbacks, the range in pCO 2 between the strongest run-off feedbacks and no feedback being 32 ppm 100 kyr after a 5000 GtC emission, compared with 81 ppm for temperature feedbacks.After 1 Myr, the model runs with strongest run-off-weathering feedback are only back to 291 and 345 ppm respectively for 1000 and 5000 GtC emissions (compared with 275 and 288 ppm respectively for the strongest temperature-weathering feedback).
The effect of whether run-off-dependence is explicit (Eqs.10 and 11) or implicit (Eqs.10, 11 and 12) is relatively large.When both carbonate and silicate weatheringrun-off feedbacks are switched on, there is a 44 ppm difference between the two in atmospheric pCO 2 1 Myr after 5000 GtC of emissions.This is largely down to differences in k run ; k run ≈ 0.02 for GENIE output, whereas k run = 0.045 was used for the implicit calculations, following Berner and Kothavala (2001).Note that the GENIE values of k run are closer to the real world global average value of k run ≈ 0.02-0.03(New et al., 1999;Fekete et al., 2000;Fekete et al., 2002).However, both the GENIE and the data-derived values of k run come with the caveat of very low linear correlations between temperature and run-off (see Sect. 3.4 of Colbourn, 2011).
Ensembles with both run-off-weathering feedbacks (R_Ca and R_Si) on, covering the published ranges for the constants β and k run are shown in Figs.11 and 12. Like with the overall effect of having the run-off feedbacks switched on or not, ranges in atmospheric pCO 2 values between end members are modest after 100 kyr even in the 5000 GtC scenario.The difference between β = 0.48 and β = 1.12 is 16 ppm; the difference between k run = 0. these ranges to be increasing through to year 500 k, where they become 43 ppm and 69 ppm, respectively.

Productivity dependence of weathering
Following Lenton and Britton (2006), weathering fluxes are given an optional (arbitrary) linear dependence on land plant/biosphere productivity (P): Productivity can be measured as either net or gross primary productivity (NPP or GPP), taken explicitly from a coupled land vegetation scheme such as ENTS (Williamson et al., 2006).However, difficulties were experienced in obtaining a reasonably calibrated model setup using ENTS in GENIE, hence an option to parameterize productivity as a function of pCO 2 following the approach of Berner (Berner, 1991) was explored: where C is atmospheric pCO 2 .This is a widely used approximation of the fertilisation effect of CO 2 on land plants.
The effect of the productivity feedback is shown in Fig. 13.As with the run-off feedback, there is a significant effect only when the emissions pulse is large (5000 GtC).Again, the effect is significantly smaller than that due to temperature feedbacks, the range in pCO 2 between the strongest productivity feedback and no feedback being 47 ppm 100 kyr after a 5000 GtC emission, compared with 81 ppm for temperature feedbacks.

River routing
Weathering fluxes are routed to the coastal ocean using simple continental "roof" routing based on the average corresponding real-world altitude of each land grid cell, and simulated topological network (STN) data (Vörösmarty et al., 2000b, a).There are two river routing schemes available.
The first is simple roof routing (using the average altitude of each land grid cell).Here, land grid cells are marked with numbers corresponding to the four cardinal directions (N, S, E, W), and a path is followed from each land cell to a coastal ocean cell.
The second uses the detailed STN routing data, which is similarly marked, but with the addition of the intercardinal directions (NE, NW, SW, SE).Roof data is used for run-off that ends in grid cells that are designated as land on the GE-NIE grid; this is about half, as due to the low resolution of the GENIE grid -36 × 36 cells for the whole globe -coastal land cells can typically cover areas that are half ocean in the much higher resolution grid of the STN data -0.5 • × 0.5 • .The higher resolution of the STN data means that for each GENIE land grid cell, water is routed to a number of different coastal ocean cells.A table containing a line for each GE-NIE land cell was produced; each line lists destination cells and the corresponding fraction of the starting land cell's water that goes to each of them.This table is read in by the RokGeM model, and used to route weathering fluxes from land to ocean.For scheme 2, coastal ocean run-off flux is scaled relative to the total land area, so that run-off flux is conserved.This factors in endorheic run-off (flux that drains  to inland seas and doesn't make it to the global ocean; mainly in central Asia, Africa and Australia -some 18 % of the land surface is covered by endorheic basins (Vörösmarty et al., 2000b, a)).
Figure 14 shows the distribution of coastal ocean flux for both routing schemes assuming constant run-off across the land surface.Scheme 2 -being partly based on the STN data -shows the most realistic distribution.One can clearly see the mouths of some of the world's major rivers, including the Mississippi, the Amazon, the Congo, and the Nile.
Note that Antarctica isn't included in the STN run-off routing data, whereas it is in the simple roof routing.At first glance it would seem sensible not to include Antarctica as there is no significant river drainage contributing to fluxes of weathered minerals reaching the ocean.However, there is still weathering happening on Antarctica, albeit bacteriamediated and sub-glacial (Sharp et al., 1999).This only applies to the simple global average weathering rather than the explicit 2-D schemes (see below) as there is no weathering contribution from ice covered areas in the latter.
When comparing the different river-routing schemes (Fig. 15) one notices little difference to results using the two schemes.It should be noted that the timescales over which the global ocean mixes are relatively short when compared to those over which the weathering-climate feedbacks operate.This negates the need for high fidelity in the river routing of weathered mineral fluxes.

Run-off dependence of weathering (2-D schemes)
For the 2-D weathering schemes, land lithological data (Amiotte- Suchet et al., 2003;Gibbs et al., 1999) from a 1 • resolution map was placed onto the lower resolution GENIE grid, this meant that some land ended up in the GENIE ocean.This was ignored (resultant weathering fluxes were "truncated" to the GENIE land grid) on account of it not being possible to get accurate values of T, R and P (which weathering fluxes depend on) in the ocean grid cells of the model.Likewise, some water ended up on the GENIE land grid.The water fractions of land grid cells do not contribute to weath-ering flux.Weathering fluxes are calibrated by scaling to the total number of GENIE land cells.

The Global Erosion Model for CO 2 fluxes (GEM-CO2)
The Global Erosion Model for CO 2 fluxes (GEM-CO2) (Amiotte Suchet and Probst, 1995) (1999); note that shale and sands and sandstones contain significant amounts of carbonate rock.The global lithological map was placed onto the GENIE-1 grid of 36×36 equal area longitude-sine(latitude) cells.Each grid cell was assigned a value for each rock type (and also water and ice) corresponding to the amount of that rock type present in it.The values for all rock types, water and ice were normalised to sum to 1 for each grid cell.Figure 17 shows the re-gridding.Note that due to imperfect grid-matching between the data sources and the model, some water appears on the GENIE land grid (Fig. 17c, l); this land does not contribute to weathering flux.Likewise land does appear in the GENIE ocean grid (Fig. 17t, similar for GEM-CO2 but 5 %); this land is truncated and not used.After truncating data to the GENIE grid, values are calibrated by scaling to the total number of GENIE land cells.

The Gibbs and Kump Weathering Model (GKWM)
The Gibbs and Kump Weathering Model (GKWM) (Bluth and Kump, 1994;Gibbs and Kump, 1994;Gibbs et al., 1999) is similar to GEM-CO2, but with the run-off dependence modulated by a lithology-dependent fractional power β as well as a linear coefficient k (Eq.15).There are also only five rock types; for the purposes of using the same lithological data (Amiotte- Suchet et al., 2003), the granite class is counted as the sum of the shield and acid-volcanic classes from GEM-CO2.A significant fraction of the land surface is designated "complex"; this signifies none of the main rock types being dominant; for the purposes of the weathering model here described, the complex class was divided into the others in the following fractions: carb=0.172, shale=0.330 sand=0.158, basalt=0.079, acid=0.017, shield=0.243, gran-ite=0.260 (Cirbus Sloan et al., 1997).Gibbs and Kump also produced a global lithological dataset, which can be used in lieu of the GEM-CO2 dataset in RokGeM (see Fig. 17b, l-t). www

Visualising the functions
The run-off dependencies of weathering flux are described by Eq. ( 15) and Table 2 for each of the 2-D schemes, where F Ca 2+ is the weathered flux of Ca 2+ (Mmol km −2 yr −1 ), R is run-off (mm yr −1 ), and k and β are dimensionless constants; they are plotted in Fig. 16.
Note that the factor of 1 2 comes from the fact that in the original papers the weathering flux is described as bicarbonate, whereas for the purposes of this study it is measured as Ca 2+ .The functions (with the Globavg scheme for comparison) are illustrated along with the temperature feedback of Eqs. ( 6) and (8) in Fig. 18.Despite their differing formulations, both of the 2-D weathering schemes produce similar results.
Figure 19 shows the distribution of terrestrial weathering flux over the land and the ocean.Areas of highest flux are the eastern side of the Americas in the tropical regions; West, East and South-East Africa; and South, East, and South-East Asia.Of the small differences that appear between the two different schemes (of order 10 %), the greatest are in the same regions of high flux.

Sensitivity tests of 0-D and 2-D schemes
Comparing the spatially explicit (2-D) weathering schemes GKWM and GEM-CO2 with Globavg (0-D) for the different weathering feedbacks, there is a noticeable difference when f_Si is switched on (Fig. 20).Both of the 2-D weathering schemes show similar results; 100 kyr after a 5000 GtC emissions pulse there is a 4 ppm difference in atmospheric pCO 2 between the two, whereas there is a 20 ppm difference between their mean and 0-D Globavg. in this instance the 0-D model starts with prescribed carbonate and silicate weathering fluxes that match the average of the 2-D models, but when matched to sediment dissolution rates during spin-up (Sect.2.4), these fluxes diminish to this lower value).This disparity in weathering fluxes between the 0-D and 2-D schemes is a consequence of the prescribed 0-D weathering not factoring in shallow-water carbonate weathering, whereas the 2-D models do.Ideally, a better comparison would be between prescribed 10 Tmol yr − 1 0-D weathering and suitably scaled down 2-D weathering.
The largest difference between the 0-D and 2-D schemes in the 5000 GtC emissions scenario occurs at year 5800 and is 113 ppm.The difference is largest around this period because for the 0-D scheme, the CaCO 3 sediments almost completely dissolve, whereas for the 2-D schemes, due to a higher preindustrial starting point, sediment dissolution does not become saturated.For the 1000 GtC scenario, where the sediments do not come near to complete dissolution for either the 0-D or 2-D schemes, there is no pronounced peak in discrepancy between the two.
Runs not including f_Si show a smaller difference (Fig. 21).Although the f_Ca-only runs show that f_Ca continues to act for the GKWM throughout the simulation, the Globavg run levels off after ∼ 50 kyr.
With 2-D weathering, the ocean carbonate system reaches different equilibria at spin-up to those of the 0-D weathering model.Surface ocean pH varies by up to 0.018 units between the 0-D and 2-D schemes (Globavg = 8.149; GKWM = 8.163; GEM-CO2 = 8.167).

Lithology dependence of weathering
As a test of sensitivity of the spatially explicit weathering in the model, different lithologies were tried using the GEM-CO2 weathering scheme.Data compiled by the teams behind GKWM (Gibbs et al., 1999) and GEM-CO2 (Amiotte- Suchet et al., 2003) were compared, along with their averages (where each grid cell contained an amount of each rock type equal to its global average) and monolithic scenarios where the entire land surface was set to a single rock type.
Results are shown in Fig. 22; for the individual lithologies, colours follow those in Fig. 17.There is little difference between the datasets of GKWM (the default; dark green) and GEM-CO2 (light green line overlapping dark green line).(c-t).Rock types are colour coded as in (Amiotte- Suchet et al., 2003).For the GENIE grid, each rock type fills a certain fraction of each cell; increasing fractions from 0 to 1 are shown by fading in colour from white (0) to full colour (1).See text for definition of "complex".
There is also little difference between the GKWM lithology and its average (dark green line vs. dotted dark green line); this is perhaps down to the influence of the "complex" class providing an averaging effect.However, there is a significant difference between GEM-CO2 and its average (light green line vs. dotted light green line).The fact that the average lithologies give higher pCO 2 values indicates that the more weatherable rock classes (carbonate, shale, basalt) are concentrated in the zones of higher weathering potential (higher temperature and run-off).
There are very large differences between the monolithic end members; shale having the greatest effect, and sandstone the least on global weathering fluxes and hence drawdown of atmospheric pCO 2 .The monolithic lithologies closest to reproducing the results of the full model are acid volcanic rock and granite, which closely follow the GKWM and GEM-CO2 full lithologies (the magenta and purple lines of the acid and granite members are in fact hidden underneath the green lines of the GKWM and GEM-CO2 members).In the GEM-CO2 scheme used, granite is not a distinct lithology, and so it is made up of half acid volcanic rock, and half shield rock.The red line shows that a monolithic shield rock world weathers slower, so the granite world must be dominated by its acid volcanic rock component.The carbonateonly world (mono carbonate) weathers quickly to begin with, but then slows after the "terrestrial neutralisation" carbonate weathering-feedback phase.The full model (GKWM or GEM-CO2 lithologies; solid green lines) overtakes the carbonate world (blue line) after year 14 k, when the silicate weathering phase is dominant.This is on account of the carbonate end member containing only a small fraction (7 %) of silicate rock.

Summary of RokGeM model description
RokGeM calculates weathering fluxes of alkalinity and dissolved inorganic carbon (DIC) dependent on, and in feedback with, inputs of land temperature (T ), run-off (R) and productivity (P ).Carbonate and silicate weathering fluxes of calcium ions (F CaCO 3 and F CaSiO 3 respectively) take the form: where 0 denotes an initial value in the feedback ({T , R, P } = {T 0 , R 0 , P 0 } for switched off feedbacks), k Ca is a constant, E a the activation energy of the silicate weathering reaction, and β (0 < β < 1) is a fractional power dependent on lithology.Productivity was parameterized following Eq.( 14).Fluxes are worked out as a global average for the basic 0-D implementation of the model, and individually for each grid cell, for the spatially explicit 2-D version of the model.In the latter 2-D case, each grid cell is apportioned between 5 or 6 distinct lithology classes.Fluxes are routed to the coastal ocean, where they interact with ocean, atmosphere and sediment biogeochemical cycles.For the temperature dependence of silicate weathering we use an activation energy of E a = 63 kJ mol −1 , based on laboratory experiments (Brady, 1991), but also explore a range of E a = 45-103 kJ mol −1 , based on field data (West et al., 2005).Recent studies based on field measurements suggest low values of E a = 49 kJ mol −1 for granites (Oliva et al., 2003) and E a = 40 kJ mol −1 for basalts (Dessert et al., 2003), at and below the bottom end of the range from West et al. (2005).These lower activation energies will produce a smaller increase in silicate weathering flux for a given increase in temperature (Fig. 5), causing atmospheric pCO 2 to recover more slowly for a given anthropogenic perturbation, although the effect is fairly modest (Fig. 7).The assumed increase of carbonate weathering with temperature (Eq.6), based on field measurements, is surprising because carbonate solubility decreases with temperature.Hence it may be capturing coeval changes in run-off or vegetation.An alternative approach would be to calculate carbonate weathering assuming run-off is equilibrated with calcite, using plant productivity to estimate soil CO 2 and using air temperature to calculate solubility.However, the model results are not very sensitive to the existence or not of a temperature dependence of carbonate weathering (Fig. 6).

Geosci
The optional dependencies of silicate and carbonate weathering on CO 2 , following a Michaelis-Menton dependence of plant productivity on CO 2 (Eq.14), neglects the sen-sitivity of plants to the water cycle.For example, rising CO 2 tends to trigger stomatal closure, increasing plant water use efficiency, which may reduce evapotranspiration and increase run-off (Beaulieu et al., 2010).However, even when our relatively crude and strong sensitivity of weathering fluxes to CO 2 is switched on, its effects on atmospheric pCO 2 are smaller than those of temperature feedbacks (Fig. 13).

Additional factors affecting weathering
With weathering, there are many variables to consider, which have a high degree of interdependence.As well as temperature and run-off, there are also the biotic components vegetation and soil to consider.These are to some degree summarised in the basic productivity dependence outlined in Sect.3.4, although there is the risk of double-counting in that the temperature and run-off dependencies arise from field data on an Earth that is clearly biotic, and so presumably already include the effects of biology.Further nuances of weathering lead to the consideration that perhaps chemical and physical weathering processes should be separated out.Erosion has been shown to be a major influence on weathering (West et al., 2005) and differing levels of erosion divide weathering into "transport-limited" and "kinetically limited" regimes, which should also perhaps be separated in any weathering model that strives for realism.Physical weathering is also highly dependent on wind (dust abrasion) and relief, which is in turn a function of altitude.Wind driven dust is also a source of minerals that add to riverine flux usually deigned to be from weathering (Hilley and Porder, 2008).Another factor is ground frost, which leads to cracking and, therefore, mineral release.
A recent study of weathering in Japan (Hartmann et al., 2009) emphasises the role of slope as a factor influencing weathering flux.Weathering is inversely correlated with slope; with increasing slope, reaction rates per volume of water with mineral surfaces decreases.
Also in need of consideration are human factors.Features of the Anthropocene (Zalasiewicz et al., 2008) particularly pertinent in their effects on carbonate and silicate weathering processes include large increases in erosion rates due to unsustainable agricultural practices; the production of acid rain; the altering of vegetation types and areal cover; and the addition of mineral fertilisers to soil.These may be specified as boundary conditions, but applying them in a timedependent manner over model runs spanning decamillennia (10 4 ) to lakhs (10 5 ) of years would involve pure speculation, considering global macro-economic effects cannot be predicted even years (cf. the "credit crunch") or centuries (cf. the industrial revolution) in advance.
Unfortunately, many of these factors will likely never be adequately parameterized in models on the scale of today's EMICs.This is on account of their variability over small spatial scales.Even across single kilometer-scale watersheds there can be wildly different weathering regimes; mountainous regions where erosion is dominant, merging into forested regions downstream, where thick soils mediate the weathering reaction.
Progress has been made recently in formulating general weathering rate relationships (Lasaga and Lüttge, 2001;Lüttge, 2006;Hellmann et al., 2010).Through looking at the molecular surface interactions involved, and using equations analogous to those for crystal growth, general relationships based on Gibbs Free Energy and separate near-and far-fromequilibrium regimes were determined.Attempting to parameterize these for use in an EMIC proves to be problematic, however, on account of variables relating to concentrations in land water of the necessary aqueous carbonate species not being readily modelled.
A mechanistic model of weathering has recently been developed however (Roelandt et al., 2010).The Biosphere-Weathering at the Catchment Scale (B-WITCH) models soil geochemistry and upscales it to a 0.5 • × 0.5 • continental scale.Chemistry is modelled directly through kinetic laws derived from Transition State Theory.Carbon and water inputs are provided through a coupling to a dynamical vegetation model.This mechanistic approach has the advantage over the parametric approach employed in RokGeM (and the other global models it was based on) in that it captures more detail on both temporal and spatial scales.Global parametric weathering models can only capture broad features relating to climate (temperature and run-off).Mechanistic models such as B-WITCH however, are more able to capture the rapid changes of the Anthropocene (Beaulieu et al., 2012), by modelling smaller scale changes in soil composition and chemistry.However, this comes with a computational overhead, and so far B-WITCH has only been used to model subcontinental scales.
Ideally, a number of things would need to be in place within the GENIE model for a weathering module to provide realistic results.Accurate temperature and moisture fields for the continental interiors would be important, as would wind fields and orography for physical weathering.If it is the case that a 3-D atmosphere (i.e. the IGCM) is necessary for this, then the feasibility of very long runs (10 5 -10 6 yr) becomes an issue as the IGCM 3-D atmosphere module runs nearly 100 times slower than GENIE-1.

Errors in modelling
As documented in Colbourn (2011) Sect.2.5, a significant source of error in the RokGeM model is the quality of the model fields produced by the wider GENIE model (in particular the EMBM) used as inputs.The ideal solution of course would be for these inputs to be better modelled (or simply better tuned) so as to get more accurate (in comparison to real-world data) results.Short of this -which is perhaps an unrealistic expectation given the very coarse resolution of the model as configured -an option to calibrate input fields to data could be created and used.
Another source of error is the run-off-as-a-function-oftemperature parameterization (see Sect. 3.3,Eq. 12).Using relations from box-modelling (Berner and Kothavala, 2001), the option is available in RokGeM to use temperature as a proxy for run-off.This might be considered a windfall, in light of the fact that run-off fields are perhaps even worse (less accurate) than temperature fields in GENIE (as configured; see Sect.2.5 of Colbourn, 2011).However, the correlation between run-off and temperature in real world data seems tenuous enough to stretch the credibility of the parameterization.It is interesting to note that in the archetype WHAK weathering model (Walker et al., 1981) Hartmann et al. (2009) discusses the role of dissolved silica (DSi) -a product of silicate weathering -in providing nutrients for ocean-dwelling diatoms, part of the "biological pump involved in the sequestration of atmospheric pCO 2 .High tectonic activity in the Pacific "ring of fire" leads to high levels of weathering, a high DSi flux entering the ocean, and hence significant CO 2 sequestration.The ocean biogeochemistry module of GENIE, BioGeM, includes a silica cycle already, so this should be a relatively straight forward addition to RokGeM.

Higher resolutions
The GENIE model is capable of being run at a number of different spatial (and temporal) resolutions and grids (Marsh et al., 2011;Lenton et al., 2006Lenton et al., , 2009)).The option to run RokGeM at different resolutions is included by way of a parameter which specifies the directory that contains the lithological data files used.The GKWM and GEM-CO2 lithological data (see Sect. 3.6.1 and Sect. 3.6.2) were also gridded onto a 72 × 72 equal-area longitude-sine(latitude) grid (double the default resolution), but this has not been tested on account of other climate and biogeochemical model configurations that are yet to be done, and the longer run times involved in running the model (it already takes ∼ 5 weeks of CPU time to run 1 Myr with the 36 × 36 version; the 72 × 72 version would run 8 or even 16 times slower on account of it also including the 16 level ocean).
A version of GENIE using a 16 level, rather than 8 level (GOLDSTEIN) ocean is now being widely used (Ridgwell et al., 2007;Goodwin and Ridgwell, 2010) on account of its higher-resolution ocean better reproducing biogeochemical tracer distributions.This version, which is only twice as slow to run as the version used for all results previously presented, is currently undergoing testing using RokGeM.With the current latest CPUs, it is now likely feasible to run a full carbon-cycle version of GENIE (including RokGeM, which only takes ∼ 3 % of the model run time) at 36 × 36 × 16 resolution, for 1 Myr, inside a month.

Applying the model to glacial-interglacial cycles
There is potential for the anthropogenic perturbation of the carbon cycle to cancel the next ice age (Archer, 2005), or even flip the Earth System into an ice-free "Hothouse" climate, last seen in the far-distant past, 3Ma in the Pliocence.This poses the unresolved question of just how large the perturbation has to be for this to occur.Mysak (2008) reviews glacial inceptions.He concludes that due to Milankovitch forcing alone, the next glacial inception is unlikely to be within 50 kyr, as long as atmospheric pCO 2 remains above preindustrial levels (>280 ppm).Mysak models glacial inceptions using the McGill Paleoclimate Model (MPM), an EMIC that has detailed vegetation and ice-sheet parameterizations that allow a good representation of the important albedo feedbacks affecting the timing of glacial inception, despite its 2-D latitude-depth ocean, and compartmentalized continent grid being of lower resolution than GENIE.Other important feedbacks modeled in the MPM are the orography-temperature feedback, and freshwater fluxes and their affect on polar heat transport via the thermohaline circulation.For the case of atmospheric pCO 2 fixed at 300 ppm, it was also found with the MPM that glaciation is avoided for the totality of a 100 kyr run.Runs including a large anthropogenic warming perturbation (equivalent to thousands of GtC) show climate has little memory following global warming; i.e. the timing of far distant glaciations is unaffected.However, (Archer, 2005) argues, and results presented here suggest, that this is not the case with an interactive carbon cycle present.The long tail in the perturbation stretches to possibly subsume future glaciations.
Using the CLIMBER-2 Earth System Model, and including a parameterization of ice-sheet dynamics and a summer insolation-minimum for glacial inception, it was determined that with 5000 GtC of emissions, the next glaciation is postponed indefinitely (Archer and Ganopolski, 2005), as an atmospheric concentration of >400 ppm CO 2 persists.For 1000 GtC emissions, the next glaciation is postponed from 50 kyr to 130 kyr hence; this is in agreement with Mysak (2008), as the atmospheric pCO 2 stabilizes at a level slightly above 300 ppm in this case.The above results are without the inclusion of weathering, however.With the inclusion of the draw-down of atmospheric pCO 2 through weathering, presumably the delay in the inception of glaciation will be less.In the case of the 5000 GtC, the delay will not be permanent.That is, unless the climate reaches a new equilibrium on account of changes in the land carbon inventory.In order to resolve the question satisfactorily with (the fully interactive ESM) GENIE, in addition to the RokGeM weathering model, the key processes of ice-sheet dynamics and methane clathrate release need to be adequately parameterized.There is an option to include the 3-D ice sheet model GLIMMER (see Supplement of Lenton et al., 2007) in GENIE.The land model ENTS needs to be active in GENIE to determine whether there is a shifting in the equilibrium land carbon storage.The silicate weathering feedback timescale of ∼ 110 kyr given by GENIE-RokGeM is similar to the glacial-interglacial timescale over the past few cycles, as measured in the polar ice cores (Lüthi et al., 2008).This suggests that it could be sufficient to mitigate against severe anthropogenic disruption of the glacialinterglacial cycle.

Applying the model to Phanerozoic events
The inclusion of RokGeM into GENIE allows for long-term carbon cycle experiments to be performed in a manner general enough (geographic boundary conditions can be readily changed) to be applied to any geological interval or event.
A caveat here is that the distribution of rock types is less well known the further one goes back in Earth history.Therefore, the 2-D version of RokGeM can be used for Cenozoic and Mesozoic studies, but simulations will be restricted to the 0-D version for Paleozoic and earlier studies.Past events in Earth history involving large perturbations in climate and the carbon cycle include the Devonian rise of plants ( Le Hir et al., 2011), basalt outpourings at the End Permian (Siberian traps) and End Cretaceous (Deccan traps), and the release of fossil carbon at the Paleocene-Eocene Thermal Maximum (PETM).These events are all future possibilities for modeling with GENIE-RokGeM.
The PETM is arguably the best analogue of future anthropogenic climate change we have from the geological past (Pagani et al., 2006).Previous work with GENIE (Panchuk et al., 2008) used a sediment data-model comparison, and modeled ocean circulation patterns resultant from relevant paleogeography, to constrain the source and amount of carbon released during the PETM.A pulse of 6800 GtC into the atmosphere (and/or ocean) is modeled as the minimum required to match the sedimentary records.The version of GENIE used by Panchuk et al. (2008) has a fixed weathering flux applied uniformly throughout the global ocean.This is similar to the 0-D version of RokGeM described herein, but without the coastal placement of weathering flux.Switching to coastal input only minimally slows down oceanic uptake of atmospheric carbon.However, the switching on of active feedbacks, and 2-D weathering, has a large impact, even on the 10 2 -10 3 yr timescale relevant to initial "upward" perturbation of the carbon system following the PETM (see Fig. 21).This suggests a limit somewhat lower than 6800 GtC as a lower bound for atmosphere-ocean input of carbon at the PETM.

Conclusions
A new spatially explicit weathering model, RokGeM, has been developed, designed as a modular component of a wider Earth System Model (GENIE).Following earlier modelling (GKWM, Bluth and Kump, 1994; GEM-CO2, Amiotte- Suchet et al., 2003) it has dependencies on run-off and temperature.It also has an optional explicit dependency on biological productivity, which future work will link to the land surface component of GENIE, ENTS (Williamson et al., 2006).Other works in progress include the better optimising of input fields (of temperature, run-off and productivity), and the extension of the model into higher-resolution and palaeo configurations.Not to mention the many factors complicating the accurate modelling of carbonate and silicate weathering on a global scale.Nevertheless, RokGeM is the first spatially explicit weathering model to be fully dynamically interactive with the carbon and climate cycles and able to integrate out to steady state over a million years.svn co https://svn.ggy.bris.ac.uk/ subversion/genie/branches/greg-s_branch? r6324 -username=genie-user Careful when copying and pasting to make sure all characters come out, and there are no spaces in the URL.You will now have a copy of the model in a directory named greg-s_branch?r6324.
Rename this directory genie.
2. You will now need to edit the files genie-main/user.shand genie-main/user.mak in order to set environment variables (directory paths, compiler options and library locations) relevant to your system.Possible options you may need are commented out with a #.Replace ∼ /genie_dev and $(HOME)/genie_dev with the path to the directory containing the model (directory genie from step 1 above).

Test the compilation of the model by typ-
ing make from the directory genie-main.
Most likely it won't work first try and will need some tweaking of variables and flags.See https://source.ggy.bris.ac.uk/ wiki/GENIE:Compiling (username: genie-user, password: tosca) for tips on compilation.5. Assuming all tests pass, move on to the following section to learn how to run experiments as done for this paper.

A1.3 Automating splitting long runs with multiple stages into manageable pieces
1. Copy the contents of genie-tools/genie_rootdirstuff into your root directory (the directory where the genie directory resides).This contains the following: -genie_output.
-genie_forcings (this contains all the forcing files necessary for the experiments shown in this paper, and more).
To be useable as individual runs (to run using the method in A1.3 below), these need to be unpacked using Mathematica (see A1.6).
-genie_runlog (to store the terminal output from script runs).
-genie_archive (to store zipped archives of results) containing directories called fresh, to_resub and to_clearout.
results (to store collated results from Mathematica output in -see below).
-$4 = NPARTS: the number of parts to the experiment -e.g. 3 for a 2-stage spin-up, then a main emissions experiment.
-$7 = MAXYEARS: your individual job length in years, e.g.4000; means the job is broken into chunks of 4000 yr each (this takes < 5 h on the UEA ESCluster).
-$9 = MINJOBTIME: the crash tolerance -min number of seconds allowed between resubmits before job is killed (10 is good for testing, 60 for actual runs).Hopefully, on screen in the Mathematica notebook will be displayed the bottom of the most recent run logs -which should give you an idea of why runs have crashed, if any have.If runs have crashed, a script will be generated in genie_runlog/expts called deadjobs_XXXXXXXXXX (where the Xs represent a date string); run this script to automatically resubmit crashed runs.Note that if you want to start runs from completely afresh, it is best to delete their directories first in genie_output and genie_archive.Leaving these in place can interfere with the processing of this part of the script.
5. Once runs are finished, or on their way, execute "Collate Results:" to import data from time series files and export to .xlsfiles (i.e.collate the ensemble time series into spreadsheets).Results (and any graphs created by setting saveallgraphs=1) will be saved in folders in the results folder corresponding to ensemble member names, experiment parts and date produced (if sameoutdir = 0 -see multi_ensemble_params.nb).Options (divisions, saveallgraphs, timescaleanalysis, timescalefitting, savenetcdfmovies -see multi_ensemble_params.nb)allow different combinations of runs to be collated, graphs to be plotted and saved, timescale analysis to be performed (specific to emissions runs where CO 2 is falling), tables to be created and exported along with graphs to LaTeX form, and animated gifs of various forms to be produced from netcdf output.6. Execute sections "a.Set-Up" and "b.Interactive time series plotting" in the section "Interactive Plotting:" to get an interactive time series plotting widget with a fairly self-explanatory GUI.The tables displayed below the graph will be of use in picking which lines (ensemble members) to display (tab).To save a graph, check and then uncheck the "save graph" box (leaving it checked will produce a new file every minute)."c.Interactive timescale analysis plotting" is similar, but with timescale analysis (e-folding curves) for global average pCO 2 and surface temperature.
7. For interactive netcdf output, see section "d.Net-CDF".Execute section "i.Set-up" and then have a play with the other sub-sections, which get more complicated (and computer resource intensive) as they advance.Each one produces a GUI with geographic output that should be fairly intuitive.Note that if there are more (or less) than 2 ensemble variables, then the code needs to be altered where the grey commented out bits are.For example:

Fig. 1 .
Fig. 1.Schematic of the relationship between the different model components comprising the GENIE-1 model.Outlined in bold is the weathering module RoKGeM described in this paper, while the ocean-atmosphere-sediment carbon cycle and climate model (Ridgwell and Hargreaves, 2007) is delineated by a dashed box. Figure adapted from Ridgwell and Hargreaves (2007).
Fig.2.Sedimentary carbonates moving toward equilibrium in model spin-ups.All ensemble runs are shown, covering the full range of parameters explored in this paper.Note that most of the spread (<30 & >35 wt% CaCO 3 ) is from ensemble members with spatially explicit weathering and greatly differing lithologies (see section §3.7); the observed real-world value is 34.8 wt%(Ridgwell, 2007).After ∼50 kyr, sedimentary CaCO 3 by wt% has levelled off for most runs (a).Some are slower than others to reach equilibrium.The vast majority of runs have CaCO 3 wt%s less than 1% different to their final values at 100 kyr after 50 kyr (b; there are 2 outliers not shown).After 80 kyr, all runs runs are less than 1% from their 100 kyr values (d).By 90 kyr, all but one run has a wt % of sedimentary CaCO 3 less than 0.1% away from its value at 100 kyr.At this stage, sedimentary CaCO 3 is settled at equilibrium, only fluctuating a few parts per ten-thousand over decamillenia.

Fig. 3 .Fig. 3 .
Fig.3.Time-series visualisation using an example result of atmospheric pCO 2 (dependent variable) over 1 Myr for two different emissions pulses: linear (a, b, c) and logarithmic (d, e, f) scales for dependent variable with a linear timescale (a, d), a logarithmic timescale (b, e), and multiple linear timescales (c, f); multiple linear scales for both the dependent variable and time (g); separate plots for each emissions scenario (h) [continues on next two pages].

Fig. 4 .Fig. 4 .
Fig.4.Evolution of atmospheric CO 2 over 1 Myr for 1000 GtC (bottom) and 5000 GtC (top) emissions pulses, vs. a control run with no emissions (grey lines at bottom of plots), with short-circuiting of the atmosphere on/off.

Fig. 12 .Fig. 12 .
Fig. 12. Evolution of atmospheric CO 2 over 1 Myr for 1000 GtC (bottom) and 5000 GtC (top) emissions pulses, vs. a control run with no emissions (grey lines), with different values of runoff-temperature correlation constant.

Fig. 13 .Fig. 14 .
Fig. 13.Evolution of atmospheric CO 2 over 1 Myr for 1000 GtC (bottom) and 5000 GtC (top) emissions pulses, vs. a control run with no emissions (grey lines), with carbonate and silicate weatheringproductivity feedbacks (P Ca and P Si respectively) on/off.

Fig. 17 .Fig. 17 .
Fig.17.Distribution of lithologies from GEM-CO2 (Amiotte-Suchet et al., 2003) (a, c-k)  and GKWM(Gibbs and Kump, 1994) (b, l-t); original source data is shown in the large panels (a, b; GEM-CO2 at 1 • resolution and GKWM at 2 • resolution); data interpolated onto the GENIE grid (36x36 equal area longitude-sine(latitude) cells) is shown in the smaller panels (c-t).Rock types are colour coded as in (Amiotte-Suchet et al., 2003).For the GENIE grid, each rock type fills a certain fraction of each cell; increasing fractions from 0 to 1 are shown by fading in colour from white (0) to full colour (1).See text for definition of "complex".

Fig. 19 .Fig. 19 .Fig. 20 .Fig. 20 .
Fig. 19.Distribution of weathering flux on land (grey-scale), and once routed to the coast, in the ocean (colour-scale) for a post-spin-up pre-industrial climate, for each of the 2D weathering schemes, GEM-CO2 (a) and GKWM (b).The GEM-CO2 lithological dataset is used.There is a small difference between the schemes, as shown in c (land weathering rate) and d (coastal input rate).

-
Fortran and C++ compilers -svn version control software (necessary to download the model) -netCDF for producing model output (necessary for successful model compilation) -Grid Engine software (i.e.be on a cluster) -Mathematica (for ensemble generation and data processing) A1.2 Obtaining the model and getting it functioning 1.In the terminal, execute the following line to download the model from the svn repository:

4.
Once the model compiles without error, run the tests to make sure it is functioning properly.Again from the directory genie-main, type: make assumedgood [this makes the assumed good files that the test results are compared to] make test [runs a basic test] make testebgogs [tests the atmospheric and ocean physics parts of the model] make testbiogem [tests the biogeochemistry module].
config are stored.Again, config files relevant to the experiments shown in this work are contained here.The configs for the runs shown in the Figures are as follows (Figure, ensemble_0*.csv):

Table 2 .
Constants for 2-D weathering functions following the formula in Eq. (15) with F in units of Mol and R in units of m yr −1 .Note: shield and acid lithologies are combined into Granite in GKWM.
constructing a spatially explicit model.The types of rock, with their rates of weathering, are shown in Table2.The fractions of each rock type given to weather as either carbonate or silicate rocks is shown, followingGibbs et al.