Insights into biogeochemical cycling from a soil evolution model and long-term chronosequences

Despite the importance of soil processes for global biogeochemical cycles, our capability for predicting soil evolution over geological timescales is poorly constrained. We attempt to probe our understanding and predictive capability of this evolutionary process by developing a mechanistic soil evolution model, based on an existing model framework, and comparing the predictions with observations from soil chronosequences in Hawaii. Our soil evolution model includes the major processes of pedogenesis: mineral weathering, percolation of rainfall, leaching of solutes, surface erosion, bioturbation, the effects of vegetation in terms of organic matter input and nutrient cycling and can be applied to various bedrock compositions and climates. The specific properties the model simulates over timescales of tens to hundreds of thousand years are, soil depth, vertical profiles of elemental composition, soil solution pH and organic carbon distribution. We demonstrate with this model the significant role that vegetation plays in accelerating the rate of weathering and hence soil profile development. Comparisons with soils that have developed on Hawaiian basalts reveal a remarkably good agreement with Na, Ca and Mg profiles suggesting that the model captures well the key components of soil formation. Nevertheless, differences between modelled and observed K and P are substantial. The fact that these are important plant nutrients suggests that a process likely missing from our model is the active role of vegetation in selectively acquiring nutrients. This study therefore indirectly indicates the valuable role that vegetation can play in accelerating the weathering and thus release of these globally important nutrients into the biosphere.


Introduction
Soils play a major role in many global biogeochemical cycles due to their position at the interface between the atmosphere and lithosphere.For example, soils influence the flow of water to rivers and vegetation, they govern the flux of nutrients between the lithosphere, vegetation and rivers and are also a source and sink of gases to the atmosphere.A quantitative description of the evolution through time of the processes and properties within soils is therefore of great interest.
This study is motivated by several important global-scale questions that a dynamic soil model will enable us to investigate further.An example being the exchange of plant nutrients between the soil and vegetation, of particular importance is phosphorus which is almost completely derived from the lithosphere and considered a limiting nutrient for many tropical forests across the globe (Vitousek and Sanford Jr., 1986;Vitousek et al., 1993;Quesada et al., 2012).Another important Earth system process is the long-term carbon cycle, specifically the relationship between silicate mineral weathering and atmospheric CO 2 concentrations.Over multi-million year timescales atmospheric CO 2 concentrations are governed by the balance between silicate mineral weathering and CO 2 outgassing from volcanic and tectonic activity (Urey, 1952).Increased levels of atmospheric CO 2 promote the weathering of silicate minerals, which in turn, indirectly consumes atmospheric CO 2 (Walker et al., 1981).This weathering process which occurs within soils is also affected by many other factors such as temperature, precipitation, pH, soil depth and vegetation dynamics.The influence of each of these factors is hard to quantify from field studies Published by Copernicus Publications on behalf of the European Geosciences Union.
alone and current modelling attempts lack true process-based weathering feedbacks within soil profiles.Instead, to test theories such as the role of vegetation and silicate mineral weathering in modulating Earth's climate (Pagani et al., 2009) we need process-based models which integrate all interacting processes, including, mineral weathering, physical erosion, tectonic uplift, soil depth evolution, soil hydrology and vegetation interactions (Goddéris and Donnadieu, 2009).
Existing models of pedogenic processes are largely aimed at understanding landscape scale processes (Yoo and Mudd, 2008;Minasny andMcBratney, 2001, 1999;Dietrich et al., 1995) and largely focus on rates of soil production and pay less attention to biogeochemical processes occuring within the soils.Models developed to study whole profile evolution include those of Vanwalleghem et al. (2013); Cohen et al. (2010); Salvador-Blanes et al. (2007), which, like this study attempt to model the evolution of soil resulting from exposed bedrock over geological timescales.These models track the vertical profile of particle size distribution through time by implementing a depth dependent soil production rate, chemical and physical weathering and overturning due to bioturbation.However, these models do not include a liquid phase so chemical weathering processes or losses from the profile due to leaching are greatly simplified.
Soil models which do include such biochemical processes exist but these attempts generally focus on very specific microscale processes such as mineral dissolution and/or vegetation interactions and are not designed for understanding pedogenic processes (Goddéris et al., 2006;Wallman et al., 2005;Warfvinge and Sverdrup, 1992).An attempt to couple such processes with a pedogenesis model is the SoilGen1 model (Finke and Hutson, 2008).This model simulates the evolution of nutrient, carbon and pH profiles, however, the model requires a large number of soil properties for initialisation and can thus only predict changes in existing soil profiles.
The model which on conceptual grounds we view as having the most potential for our purposes is the pedogenesis model developed by Kirkby (1985).This model is recognised as a pioneering attempt to model biogeochemical soil processes in the context of understanding hillslope processes (Hoosbeek and Bryant, 1992;Minasny et al., 2008).The model meets the criteria of being based upon physical processes, yet is sufficiently simple to allow the mechanisms and feedbacks behind the resulting properties to be understood.
Over recent years a range of chemical and physical soil chronosequence data, a valuable means of evaluating our understanding of evolutionary processes in soil profiles, has also become available.A good example is the soils which have developed on the lava flows of Hawaii (e.g.Chadwick et al., 1999;Porder et al., 2007).However, to our knowledge, efforts to make complete use of these soil data sets and synthesise them within one consistent, process-based modelling framework have been limited.
The purpose of this paper is to introduce a soil evolution model based on the framework described in Kirkby (1985) and explore how well this updated model can reproduce current soil properties, by placing a strong emphasis on evaluation with data.Specifically we will demonstrate how the model may be used to further our understanding of longterm nutrient cycles.While the model is based on Kirkby (1985), there are many aspects of this model which need to be changed to fit our purposes.For example, to explore long-term nutrient cycles we would like to model individual nutrients, whereas the original model of Kirkby (1985) is structured such that only three profiles are simulated; a weathering profile, organic profile and inorganic profile.The components of each of these are greatly simplified, for example, solubility is expressed as a linear relationship with the proportion of original material remaining and it is suggested that pH is then also derived from this proportion.Therefore pH does not dynamically interact with solubility and is also not altered by changes in the organic profile.
In this paper we describe with equations the individual processes and mathematical basis of the updated soil evolution model.Following on from the model description, the basic performance of the model is explored.This demonstration of the model's capability is based on a hierarchy of model simulations starting with a profile subject to weathering and leaching only, with each further simulation including an additional process.We then evaluate the model with soil chronosequence data from Hawaii, demonstrating what we can learn from such a model.The focus here is on soils of tropical systems, however, the model can be applied to other biomes by adjusting the appropriate input parameters.

Model description
The process of soil evolution is conceptualised as a vertical profile of bedrock which undergoes both chemical and physical weathering resulting in an altered profile which we term soil.The formation of soil begins when water percolates into bedrock and initiates mineral dissolution.Chemical weathering in the model is based on the central assumption that dissolution equilibrium is reached between the rock minerals and percolating water (Kirkby, 1977(Kirkby, , 1985)).In the model, water acts directly on the elemental oxides of the parent material rather than on rock minerals.The oxide composition and density of the bedrock provides the initial conditions for the modelled weathering process.
The simulated percolation of rainfall through the profile provides the mechanism and rate for losses of dissolved ions from the soil layers.The modelled soil may deepen as a result of steadily increasing percolation through the profile resulting from the increasing pore space which is created by the leaching of dissolved rock oxides and also by the redistribution of soil by bioturbation and from direct removal by vegetation.The specific soil properties the model predicts are: soil depth, pore space, the proportion of the initial elemental oxides remaining in each soil layer, pH of the soil solution, organic carbon and pore CO 2 concentration (Fig. 1).The processes included in the model are, chemical weathering of bedrock elements, percolation of rainwater, leaching of weathering products, surface erosion, bioturbation, plant litter decomposition and vertical transport, CO 2 production and diffusion and nutrient cycling.As well as chemical weathering, other adopted processes from Kirkby (1985) include losses of solutes via leaching, surface erosion, biological mixing and ionic diffusion.The main way in which the model differs from that of Kirkby (1985) is that we keep the bedrock elements separate throughout model calculations.Keeping chemical elements separate allows us to explore more comprehensively the individual cycling and feedbacks of important elements.The key elements needed to understand the model rational are detailed in the following sections.

Equilibrium reactions
As already indicated, the model assumes that dissolution equilibrium is reached between the rock oxides and the percolating waters.Although a simplistic assumption, as a first order approximation, this is preferable to a formulation using kinetic dissolution equations which are particularly difficult to constrain due to the requirement of reactive mineral surface areas.Studies have also shown that the unknowns associated with kinetic reactions are very large, with weathering rates of minerals such as plagioclase behaving closer to equilibrium predictions in natural systems than to kinetic rates derived from experimental studies (White et al., 2001(White et al., , 2008)).The methods of calculating the equilibrium composition; thus the dissolution of rock oxides and subsequent pH of the soil solution are derived from Kirkby (1977) and Garrels and Christ (1965) and an example taken from Kirkby (1977) for SiO 2 is shown in the appendix.

Percolation
The rate of water flowing through each soil layer is regulated by the amount of pore space available in that layer.In the early stages of soil formation this is dependent only upon the porosity of the bedrock, however, over time, the losses due to leaching increase this porosity.The pore space is expressed as a fraction of soil volume (m 3 m −3 ) and is derived from the proportion p of original parent material remaining in the profile, where p = 1 for unweathered bedrock and p = 0 for complete loss of the original material.The total soil deficit w below depth z is calculated as and has the dimension of length (Kirkby, 1985).The coordinate system is chosen such that z is positive in the downward direction.
A simple vertical flow through the profile is assumed, with sub-surface flow resulting from the vertical variation in pore space.The percolation of water, F , at depth z is where F 0 is the rate of percolation allowed through the bedrock, and K is a site specific parameter related to hy- maximum rate of percolation, effectively occurring percolation is whichever is lowest, the maximum rate of percolation or the rate of precipitation minus the cumulative evapotranspiration from the soil surface to depth z.

Evapotranspiration
The process of evapotranspiration removes water from the soil profile.Here we calculate total actual evapotranspiration (E T ) as the minimum of potential evapotranspiration (E * T ) and mean annual precipitation (P A ): Although simple, the formulation still permits the model to operate under water stressed conditions.E * T is calculated for specific locations using a modified Hargreaves model (Hargreaves and Samani, 1985) (see appendix).This method is chosen because it requires only a small amount of climate data (temperature) for any specific location.The allocation of water loss by evapotranspiration to the different soil layers is determined by the distribution of roots through the soil profile, these are assumed to decline exponentially with increasing soil depth (Jackson et al., 1996).The e-folding rooting scale depth is z r so that the rate of evapotranspiration, E, at depth z is Rainfall minus cumulative evapotranspiration at depth z places a limit on the amount of water available for percolation: where E c (z) is the cumulative evapotranspiration from the surface down to depth z.Values for z r can be obtained from the rooting distributions compiled for different biomes by Jackson et al. (1996).When the modelled soil is shallow, the rooting depth and subsequent vertical distribution of evapotranspiration is limited by the soil depth.Rooting depth, d r , is the depth which contains a fraction, f , of the total root mass (Arora and Boer, 2003) and can be calculated by If we term rooting depth as the depth above which 95 % of the total root biomass is contained, following Arora and Boer (2003) we use f =0.9502 to aid simplicity, so that when d r is greater than the soil depth in either the early stages of soil development or in shallow soils, the above value of z r is adjusted so that d r equals soil depth.This will result in a greater proportion of roots in the top layers of soil.

Leaching
Following Kirkby (1985) the loss of solutes from the profile is calculated by mass balance, the main difference being that each individual element is treated separately.The derivation of this is provided in the appendix for completeness.

Ionic diffusion
The ions released into the solution from weathering can diffuse from regions of higher to lower concentrations, so following Kirkby (1985) we model this as a diffusive process, only again we treat each element separately.The derivation is provided in the appendix.This process is most important at the weathering front where ion concentrations are greatest and leaching losses are low (Kirkby, 1985).

Bioturbation
Bioturbation is the mixing and turnover of soil resulting from biological activity and is considered a major soil forming process (Gabet et al., 2003;Wilkinson and Humphreys, 2005;Yoo et al., 2005).Following Kirkby (1985) bioturbation is also represented in the model as a diffusive process with a depth dependent diffusivity coefficient (see Appendix).

Surface erosion
Removal of soil from the surface by mechanical processes is modelled through a denudation rate, T (m yr −1 ).Following the approach of Kirkby (1985), surface elevation, z s is lowered at a rate, dz s /dt, which is inversely proportional to the amount of original material remaining at the surface, p s , or p(z = 1).Described this way we find the following: The full derivation is provided in the appendix.Cosmogenic nuclides such as in situ 10 Be have provided measures of surface erosion for hillslope soils where soil thickness is assumed to be at steady state and thus rates of soil production from bedrock balance rates of loss due to surface erosion.Erosion rates calculated from these studies lie in the range of 10 to 100 m Myr −1 (Wilkinson and Humphreys, 2005).

Carbon fluxes, decomposition and mixing
To estimate carbon input into the soil we assume that vegetation cover is at steady state, with new carbon production equal to the losses from litterfall and root senescence.For this first presentation and evaluation of the model we simply assume a time invariant climate and annual net productivity (N P ) (kg C m −2 yr −1 ).The carbon is assigned to four different pools which are defined by the stability or turnover time of the pools.These are fine litter (e.g.leaves), coarse woody debris (e.g.branches/stems), fine roots and coarse roots and are assigned using allocation coefficients from the literature.The overall equation for the organic matter decomposition and mixing processes in the soil profile is where C i is the concentration of carbon (kg m −3 ) in pool i, the first term is the diffusive mixing of carbon through the soil profile by biological activity, k is the decay coefficient (yr −1 ) and I i is the carbon entering the soil profile from either plant litter at the surface or from root litter, which is distributed throughout the profile.The decay coefficient may remain constant with depth or decrease with increasing soil depth as observed in soil carbon studies using carbon isotopes (Veldkamp, 1994;Trumbore et al., 1995;Van Dam et al., 1997).For this study it is assumed that the decay rate, k, declines exponentially with increasing soil depth.For the fine and coarse litter, I provides a top boundary condition flux equal to α i N P where α i is the proportion of carbon production assigned to pool i.For both fine and coarse roots the input of carbon is distributed vertically throughout the profile according to Because of the much shorter timescale that these carbon dynamics operate on, we assume a steady-state carbon profile and hence solve Eq. ( 9) for ∂C/∂t = 0 and boundary conditions of ∂C ∂z = 0 at the bottom and a top boundary condition equal to the flux of carbon entering from the above litter.The carbon is not subject to the modelled surface erosion, however, given the very different timescales of the two processes this seems a reasonable simplification.A limitation of this carbon scheme is that biomass is present from the start of soil evolution rather than vegetation productivity evolving with the developing soil profile.This may result in an unrealistic vegetation enhanced acceleration of weathering in the very earliest stages of soil evolution.

CO 2 production and diffusion
Gases in soil are transported in either the pore space or in solution.Here we assume that the CO 2 produced from root respiration and from the above decomposition process is transported through the profile by gaseous diffusion only.This is modelled as a diffusion scheme: where C g is the concentration of CO 2 (kg m −3 soil air), where n is the total number of carbon pools (4 in this case)) and R c is the production of CO 2 from root respiration which is assigned from the literature and distributed throughout the profile following the same exponential function as for root carbon turnover.The effective diffusion coefficient in soil air is lower than that for bulk air due to both the smaller volumes of air filled pore space and the tortuosity introduced by soil pores.The diffusion coefficient for CO 2 is taken from Jones (1992), D c = 14.7 (mm 2 s −1 ) for 20 • C. To account for tortuosity a more realistic diffusion coefficient for soil air (D s ) is calculated using the following relationship of Penman (1940) (Hillel, 2004, p. 204).
where f a is the fraction of air-filled space, in the model this is equal to 1−p(z).0.66 represents the tortuosity coefficient, which means that the straight line path is approximately twothirds the length of the path of diffusion, so as the pore space increases the diffusive path will decrease.The CO 2 profile is also modelled at steady state so that The top boundary condition is equal to the atmospheric concentration of CO 2 and the bottom boundary condition allows no mixing out of the profile.This modelled partial pressure of CO 2 replaces the atmospheric CO 2 concentration used in the carbonate equations of the dissolution model and thus changes the charge balance of Appendix Eq. (A11), influencing the pH of the soil solution and solubilities of the rock oxides.

Nutrient cycling
Nutrient concentrations in vegetation depend on a number of factors such as the species of plant, the climate and the nutrient status of the soil.As a simplification, it is assumed in the model that the nutrients taken up, and hence those reentering the soil from plant litter and root turnover, try and meet optimum stoichiometric ratios.We use the optimum stoichiometric nutrients ratios calculated by Linder (1995) for deciduous plants.For the following elements, N : P : K : Ca : Mg : Fe, these are 100 : 10 : 35 : 2.5 : 4 : 0.2.Nutrient concentrations are calculated assuming a fixed proportion of biomass is made up of nutrients and a fixed relationship between N P and biomass production (biomass is double the mass of carbon).In the case of the soil not being able to supply enough of a nutrient to meet the optimum C : nutrient ratio, then the nutrient stoichiometry will deviate from that above.The nutrients are released into solution at the soil surface (g m −2 yr −1 ) from the fine litter pool and provide a flux surface boundary condition in the solute www.biogeosciences.net/11/6873/2014/Biogeosciences, 11, 6873-6894, 2014 transport equation (Eq.A19).The nutrients from fine root turnover are released into solution obeying the exponential decline in root distribution with depth.Although obviously not completely realistic, it has been observed that nutrients are readily lost from litter in the earlier stages of decomposition (Berg and McClaugherty, 2008) and this method makes it possible to readily incorporate the nutrients into the dissolution submodel.
Nutrients are then taken up from solution by the vegetation or leached from the system.Nutrient uptake from the soil profile is passive and controlled by the rate of evapotranspiration from each soil layer, the concentration of ions in solution in that layer and the rate of uptake required by the vegetation i.e. the fraction of biomass production calculated from N P .This process is represented by the second term in an updated form of Eq. (A19): where c n is the concentration of nutrient i in each layer (g m −3 ).Total c ni is calculated by integrating ∂E c ∂z c ni successively over each soil layer until the required annual uptake of nutrients is reached (i.e. when total uptake of nutrient i equals the nutrient production calculated from biomass production and hence turnover for the steady-state condition).R ni is the concentration of nutrient i returned from fine root turnover (g m −3 yr −1 ).When the nutrient uptake from a layer is greater than c n × t, uptake is set to c n / t.We know that plants can interact directly with soil minerals for nutrients, however, the main source of nutrients is likely the soil solution (Lucas, 2001) and this simple mechanism of nutrient uptake is employed for the first attempt at modelling longterm, plant-soil interactions.

Model solution and parameters
The model partial differential equations are solved numerically by finite-difference schemes.The leaching and denudation equations are solved by an upwind scheme (Morton and Mayers, 2005) and the diffusion equations by the semi-implicit Crank-Nicholson scheme (Morton and Mayers, 2005).The parameters used in the model runs described in Sect.3.2 are shown in Table 1.The values are selected as being the most appropriate from the literature, they are not constrained for one particular site.

Simulations demonstrating the model processes
To permit easier interpretation of our model predictions we follow a hierarchial procedure.This involves first running the model in it's most basic form and adding an additional process for each subsequent run.We can then get a clear sense of how each of the important processes influences the modelled soil properties and thus understand their importance in soil evolution within this modelling framework.
In the following simulations the oxide composition of the model bedrock is a basalt taken from the study of a Hawaiian soil chronosequence (Porder and Chadwick, 2009) (Table 2, Kona flow).For the purpose of this study the model is run first with only oxide weathering and leaching as the active processes.Other processes are then added successively in the order of surface erosion, bioturbation, organic carbon decomposition and nutrient cycling.For these simulations the profile is discretised into 10 cm deep layers and the total number of layers is chosen so that the total profile depth is greater than that reached by the weathering front during the simulation.The model time step is 0.1 year.Unless stated otherwise the model is run with the parameters in Table 1, and a mean annual temperature and precipitation of 20 • C and 1.7 m yr −1 , respectively.
The developmental state of the modelled soil profile is quantified as the proportion of each oxide remaining in the soil layers relative to that of the parent material.Values lower than one represent a relative loss from the profile compared to the initial unaltered bedrock material and values greater than one relative enrichment.A value equal to one indicates zero mobility.

Model simulations and setup for Hawaii chronosequence sites
To assess the ability of the model to reproduce real soil profiles we compare the model predictions with data from in situ soil profiles from soil chronosequences.Due to the slow nature of pedogenesis it is impossible to directly observe soil changes over these very long timescales.Instead we make use of chronosequences.These are series of soils which differ in the age of soil initiation but other factors of soil formation such as parent material, climate and topography remain constant.It is thus assumed that any differences in soil properties are related only to the differences in the age of the soil profile.
Hawaiian soil chronosequence data published by Porder et al. (2007) and Porder and Chadwick (2009), are used for comparison here.The soils have developed on volcanic lava flows on the island of Hawaii, and thus have parent material of relatively uniform composition (Table 2).Because of the wide range in eruption ages, soils from the Hawaiian island chain have been utilised in a number of studies looking at the interactions between soil age, weathering and nutrients (Vitousek et al., 1994;Vitousek and Farrington, 1997;Chadwick et al., 1999;Hedin et al., 2003;Porder et al., 2007;Porder and Chadwick, 2009).Porder et al. (2007) sampled soils on three lava flows aged 10, 170 and 350 ka, each spanning a topographic gradient and resulting rainfall gradient.Mean annual precipitation (P A ) varies from 0.57 to 2.5 m yr −1 , the highest rates of precipitation are found at the highest elevations.Mean annual temperature (T A ) increases from 16 • C at these higher and wetter elevations to 24 • C at the lower altitudes.
Consequently the sites receiving the lowest rainfall have the highest temperatures and are thus subject to the highest E T , resulting in a negative water balance (Chadwick et al., 2003).
It is important to note that the rainfall gradient has not always been this strong during the evolution of the soil profiles, this is a result of glacial periods and changes in the elevation of the trade wind inversion (Hotchkiss et al., 2000).The sites at the wet, higher elevations may have received 50 % less precipitation during most of their development, however, the temperatures during these drier glacial periods were also cooler, thus probably reducing the water lost from the profile by evapotranspiration.
The model is compared with the driest and wettest sites from each flow and an intermediate rainfall site.The following model parameters are modified from those in Table 1 to suit the Hawaiian sites.The monthly minimum and maximum and mean temperatures needed to calculate E T using the Hargreaves equation are taken from the Western Regional Climate Center (http://www.wrcc.dri.edu/).The site closest to the Kona lava flows is used and the temperatures were adjusted to T A of 16, 20 and 24 • C for the low, medium and wet rainfall sites, respectively.The estimated E * T calculated for these sites is 1.20, 1.34 and 1.48 m yr −1 , respectively.We assign productivity values for each site by relating E T to productivity using a water use efficiency (WUE) term.The WUE of a plant is the unit of carbon fixed per unit of water transpired.Assigning a WUE of 1 kg m −3 we estimate carbon productivity (N P ) values of 0.3, 0.53 and 0.48 kg m −2 yr −1 for each of the sites respectively, replicating the observed trend for Hawaiian vegetation of increasing N P with P A up to approximately 2 m yr −1 , declining for further increases in rainfall (Schuur and Matson, 2001;Austin, 2002).However, we are aware that the mechanisms behind this relationship are not the same.The decrease in the model productivity is due to decreasing evapotranspiration associated with decreasing P A , whereas, the changes in the observations are thought to be due to decreased N availability.The vertical root depth scale (z r ) is 0.26 m, the value estimated for tropical evergreen forests (Jackson et al., 1996).The bedrock oxide compositions used in the model runs are shown in Table 2.The erosion rate is set to 10 m Myr −1 because even though the soils sampled are not thought to have experienced high rates of erosion (Porder et al., 2007), even stable soils often experience erosion rates greater than 5 m Myr −1 (von Blackenburg, 2005) and values in the range of 7.7 to 12 m Myr −1 were calculated for basalts on the lip of Hawaiian volcano craters (Craig and Poreda, 1986;Kurz, 1986;Nishiizumi et al., 1990).Townsend et al. (1995) found that the turnover times of the intermediate carbon pool in Hawaii soils double with a 10 • C decrease in T A .It is unclear whether this increase in decomposition with increasing temperature follows a linear or exponential trend but here we assume a simple linear function of decomposition with mean annual temperature using the values observed by Townsend et al. (1995) to calculate the decomposition rate (k) of the coarse roots and coarse wood: The decomposition rates (k) of the fast carbon pools remain the same as in Table 1.The depth of the vertical model layers is increased to 0.25 m to improve the numerical stability of the simulations.
To determine the intensity of weathering of elements in a soil profile, element concentrations are commonly compared with those in unweathered bedrock and normalised to an immobile element such as zirconium (Zr) to give the fraction of the particular element remaining relative to bedrock (Brimhall and Dietrich, 1987) (See the appendix for a description of this method).For these Hawaiian soils Porder et al. (2007) used Niobium (Nb) as the immobile element.This provides values which can be directly compared with output from the model.
It is recognised that soils are complex systems and display a great deal of heterogeneity across even very small spatial scales.Nevertheless, we assume here that over these paedogenic timescales the soils sampled at each of these sites have been subject to the same soil-vegetation interactions.

Dissolution and leaching
With chemical weathering and leaching as the only active processes, we observe losses of the most soluble oxides in only the top 20 cm of the soil profile (Fig. 2, first column).The sequence of losses for the basic oxides is MgO > Na 2 O > CaO K 2 O (Fig. 3) after 20 thousand years of soil development.Of the non-basic oxides, we observe some depletion of P 2 O 5 and SiO 2 but very minimal losses of FeO and Al 2 O 3 (Fig. 2, first column and Fig. 3).The solubility of iron increases at lower pH values and Fig. 3  included in the model simulation and soil pH has lowered from approximately 8 to 6. Al 2 O 3 on the other hand displays reduced losses in the full simulation and this is attributed to aluminium being most soluble in either very alkaline or acidic solutions.At this stage the model displays very early signs of horizonisation, with a depleted top (or A) horizon and a slightly enriched saprolite (or B) horizon.Enrichment or deposition of the most soluble oxides at the bottom of the weathering front occurs when saturated solution from the layers above percolates into a bedrock layer with lower equilibrium solute concentrations.As discussed in Sect.5.1, this weathering sequence of basic oxides is similar to those in other studies.
The predictions based on chemical weathering and leaching processes alone demonstrate (i) an expected weathering sequence of oxide losses, (ii) very shallow soil profiles in the absence of any physical weathering or biological activity and (iii) evidence of horizonisation.

Surface erosion
The modelled process of surface erosion acts by shifting the simulated soil properties towards the soil surface, whilst removing those in the surface layers (Fig. 2, second column and Fig. 4).Erosion plays a larger role in older, more depleted soils.This is demonstrated when all processes are included in the model simulations (Fig. 4).The more weathered and depleted in original material the surface layer is, the greater the reduction in surface elevation, with a shallower profile then ensuing.Thus over long timescales surface erosion becomes an increasingly important process in the soil evolution model.If the rate of soil deepening becomes equivalent to the rate of surface denudation, soil thickness naturally reaches a steady state.

Bioturbation
Parameterised here as a diffusion process, bioturbation smoothes the oxide distributions in the surface layers (Fig. 2, third column) allowing the upward mixing of oxides from further down the profile.This action combined with that of surface erosion results in retention of mineral oxides in the surface layers (Fig. 3).Bioturbation acts to deepen the soil by removing material from deeper in the profile and mixing it into the upper layers.Bioturbation thus influences both the composition of the soil layers and the rate of soil production from bedrock.

Effect of vegetation
The addition of organic matter to the soil accelerates the weathering of all but the most insoluble oxides (Fig. 2, column 4).When carbon biomass is absent from the model simulation, soil development progresses slowly, CO 2 concentrations are equal to atmospheric concentrations, and pH remains above 7.When organic carbon is included in the model simulations pH decreases from ∼ 8 in the surface layers to ∼ 6 after 20 thousand years of soil development (Fig. 2, column 4 and Fig. 3).
This decrease in pH and increase in leaching losses is a result of the higher concentrations of soil CO 2 (Fig. 5).CO 2 concentrations increase with increasing soil depth, reaching approximately 50 times atmospheric levels where pore space is low (Fig. 5).For this simulation the flux of CO 2 out of the soil profile is within the ranges observed in studies of Hawaiian tropical forest soils (Fig. 6).
The addition of nutrient cycling into the simulation results in some retention of the oxides in the surface layers (Fig. 2 column 5), a trend also noted in the soil-nutrient studies of Jobbágy and Jackson (2001); Lucas (2001) and Porder and Chadwick (2009).However, in older or very wet soils the effect of plants on nutrient retention may be diminished due to the overriding effect of leaching losses as found by Porder and Chadwick (2009).In the model the return of basic ions in solution from litter decomposition alters the equilibrium status of the solution and slows the rate of dissolution.

Comparison of model predictions with observations
Figure 7 illustrates the performance of the model for the three different rates of annual precipitation on the young 10 ka lava flow for a selection of elements.The simulations of Ca and Na in the model are most realistic, followed by Mg, and then K and P. The model captures the slower rate of weathering losses in the driest site and higher rates in the intermediate and high rainfall sites.Mg, Ca and Na display very similar distributions of depletion in these Hawaiian soils, whereas the relative vertical distribution of Mg depletion differs from that of Ca and Na in the model.Modelled Mg weathers to deeper depths than those observed in the intermediate rainfall sites (where model N P is highest), also weathering deeper than the other model elements.Modelled K is particularly resistant to weathering compared with the observations and modelled P is even more immobile.however, the observations also exhibit little depletion of P in these young profiles.
For the 170 ka Hawi flow, both the Hawaiian and modelled soils have weathered much deeper in the intermediate and high rainfal sites compared with the younger 10 ka flow.The model captures the lower losses in the dry site relative to the wetter sites but doesn't replicate the more enriched surface layers (Fig. 9).The modelled Na and Ca profiles again match the observations most closely, reproducing the nearly totally depleted profiles at the wetter sites and even matching the depth of weathering.The depth of the Mg weathering front, however, is still too deep and the modelled K and P profiles indicate that the modelled processes are still too resistent to weathering for these elements.It should be noted that Porder et al. (2007) estimate that additions of dust to the Hawi flow averages 30 % of the total mass lost from the profiles and most of this dust is found in the top 30 cm which may obscure some of the weathering signal in these soils.
The 350 ka Pololu flow differs from the 10 and 170 ka flow by being underlain by a pahoehoe flow at 1.8 m.Pahoehoe flows are characterized by smooth, glassy surfaces and are less porous than the overlying, blocky lava flows.They therefore act as a barrier to weathering in these soils.By comparing the profiles of K with Na Porder and Chadwick (2009) show that even at this age, plants in the dry flow are still enriching the surface layers with nutrients but in the intermediate and high rainfall sites, leaching losses override any nutrient retention and the surface layers are depleted in nutrients.Figure 10 shows that for the dry site the model displays general agreement with weathering depths and again Na shows the closest match to the observations followed by Ca and Mg with K and P still too immobile at this age.The slow rate of chemical weathering of these two elements in the model also means that any depleted signal in the surface layers will also be removed by surface erosion.For the intermediate and wet rainfall sites the model captures the surface losses of Na, Ca and Mg.K is still too resistant in the intermediate site but agrees better at the wettest site.For this older soil, modelled P shows some signs of depletion but is still much more resistant than the observed profiles.The modelled profiles extend to nearly 5 m for the two wetter sites whereas the observed profiles reach a maximum of 1.8 m because of the impermeable pahoehoe layer at this depth.
Figure 11  with observations in the driest sites and for the intermediate aged, Hawi flow.Simulated pH is generally too high in the wetter sites which could be because modelled Al 2 O 3 is very insoluble (Fig. 2) so Al 3+ ions in solution may be lower.

Model processes
For the most elementary setup of the model, where weathering and leaching are the only processes included, the weathering sequence of basic oxides displays similar weathering sequences to other studies.For example Busacca and Singer (1989) observe a mobility sequence of Mg Na > Ca > K from alluvium deposits in California and White et al. (2008) observe a weathering sequence of Mg > Ca > Na > K in marine terraces also in California.For the three basic oxides found in the feldspar family of minerals, CaO, Na 2 O and K 2 O, the modelled sequence of losses follow an expected trend associated with the lower solubility of K-feldspar (or orthoclase) compared to plagioclase which incorporates the endmembers anorthite and albite (Nesbitt and Young, 1984;White et al., 2001White et al., , 2008)).Importantly, White et al. (2008) calculate that the pore waters of their chronosequence rapidly reach feldspar thermodynamic saturation and so the weathering velocity of Ca, Na and K is controlled by this thermodynamic state, the rate of which is determined by the flux of water, this being the weathering mechanism of our model.White et al. (2008) also found that the weathering of plagioclase is non-stoichiometric, i.e. there is selective removal of Ca over Na from plagioclase in their marine terraces.Thus the solubility of the oxides act independently, indicating that so far the dissolution and leaching of mineral oxides in this model is conceptually realistic.
Even though the current model does not predict secondary mineral formation, it is still possible to predict the secondary minerals likely to be present through an understanding of the sequence of minerals formed across gradients of weathering intensities.The modelled weathering sequence thus predicts the commonly predicted weathering pathway of a shift from predominantly silicate minerals such as the Mg, Ca and K feldspar family to the secondary Fe and Si containing clays such as vermiculite and montmorillonite, through to the Al and Si containing clay mineral kaolinite present in weathered soils.Eventually, in very weathered soils Al sesquioxides such as gibbsite dominate (Tardy et al., 1973).
In these poorly developed profiles, soil erosion in the model acts by effectively replenishing the surface layer with unweathered oxides from below.This action contributes to long-term biogeochemical and carbon cycles, for example, as we have demonstrated, increased rates of erosion exposes previously shielded minerals to weathering agents, enhancing chemical weathering fluxes (Millot et al., 2002;Riebe et al., 2001;Gaillardet et al., 1999).These interactions between erosion and chemical weathering have been proposed to explain periods of cooling in Earth's history (Raymo and Ruddiman, 1992;Raymo et al., 1988).To fully explore such relationships Goddéris and Donnadieu (2009) emphasis the need for a model which can track the growth of soil profiles while also being coupled to vegetation and climate.The chemical weathering and erosion that we have presented so far demonstrates that this model can provide a platform for exploring some of these theories.However, we acknowledge that erosion rates and hydrology will need to be implemented in a more mechanistic manner.Tectonic uplift, another process important over these timescales can be formulated in the model in a similar manner to the current erosion formulation.
Like erosion, bioturbation also results in greater retention of unweathered minerals in the surface layers and alters soil chemical and physical properties.We have shown that bioturbation can work to alter soil profiles at a depth greater than that reached by chemical weathering alone, greatly enhancing pedogenesis, consistent with other studies (Wilkinson and Humphreys, 2005).
The step-wise approach of including processes in the model framework has demonstrated the large influence that vegetation has on the weathering rate and subsequent development of the soil profile.The modelled vegetation affects the soil profile via three processes: by altering the vertical distribution of evapotranspiration from the soil profile, by increasing soil acidity through the production of CO 2 from root respiration and litter decomposition, and through the cycling and retention of nutrients.Of these three processes, the impact of CO 2 production on solution pH plays the largest role in our model.This work did not explore the response of soil solution pH to the pH of percolating rainwater or the increased acidity resulting from leaching losses, however, we believe that the influence of vegetation on weathering is much greater than these components.The model agrees with both field and laboratory studies exploring the effect of vascular and non-vascular plants on enhanced silicate mineral weathering rates (Moulton et al., 2000;Lenton et al., 2012).The enhanced rates of Ca and Mg weathering in the presence of vegetation supports the theory that the rise of non-vascular and then vascular plants on Earth may explain abrupt drops in atmospheric CO 2 concentrations and temperatures in the long-term records (Lenton et al., 2012;Berner, 1997).Our ability to quantify the weathering and leaching of Mg and Ca ions suggests that this model can provide a means of isolating the contribution of vegetation to weathering which is difficult to achieve in the field.

Evaluating the model with Hawaiian chronosequences
The comparison of simulated to observed elemental weathering profiles has highlighted the poor ability of the model to capture the depletion of phosphorus and potassium from the profiles.These are both important plant nutrients and required by plants in larger amounts than Ca and Mg and are thus strongly cycled (Jobbágy and Jackson, 2001).The model results for K may therefore highlight the importance of the active role of plants, mycorrhiza and faunal communities in mediating the release of this poorly mobile nutrient from minerals (Hutchens, 2009).The uptake of nutrients in the model is controlled by the rate of evapotranspiration and concentrations of the nutrient in the soil solution, however, there are a number of other mechanisms by which plants can acquire nutrients (Hinsinger et al., 2009).For example, roots can actively induce the release of non-exchangeable K from phyllosilicates by secreting H + to exchange with K.By actively taking up K from solution plants can also shift the solution equilibrium thus promoting further dissolution (Hinsinger et al., 1993;Hinsinger and Jaillard, 1993).By altering the solubility of K in our model, we show that the missing process accelerates the weathering of K from minerals by a factor of approximately 50 (Fig. 8).
The 10 ka flow is characterised by surface layers enriched in P and low amounts of P depletion in the deeper layers of the intermediate and wet sites.For the driest site (Porder and Chadwick, 2009) argue that the soils must receive additional P from exogenous sources.If the observed surface enrichment of P was due to cycling of the nutrient we would expect this surface enrichment to be balanced by depletion deeper in the profile, which is not observed.Dust can be a significant source of P to Hawaiian soils (Chadwick et al., 1999), but for these young flows Porder and Chadwick (2009) suggest that the addition of fine organic matter from nearby surroundings may explain the additions.Without these external sources of P, the relative immobility of modelled P may well be representative of these young soils.Also in the Hawaiian soils P losses are correlated with Fe losses in the old and wet sites.Fe can bind with P and may drive the losses in these lower oxygen, reducing soils (Porder and Chadwick, 2009), a process not represented in this model.
These comparisons show that for all profiles modelled plant nutrients P and K are not in agreement with the observations.That Na, which is not an essential plant nutrient, shows the best match, followed by Ca, which is a plant nutrient but is thought to be taken up in amounts equivalent to availability and depends on water flow to the vegetation (Knecht and Goransson, 2004), suggests that it is the process of nutrient uptake which the model is not reproducing realistically.The good agreement with Na, particulary in the intermediate rainfall sites where plants play an important role in nutrient distributions suggests a good model understanding of the other soil processes included in the model and thus that the model is a reliable tool for further developing our understanding of nutrient dynamics.

Limitations
We acknowledge that there are limitations and important pedogenic processes missing from the model.For example: (i) the model does not predict secondary mineral formation or size fractions so features such as cation adsorption and soil structure associated with these properties are overlooked, (ii) porosity is very simple, and pores are assumed to be free draining and connected, which for tropical soils may be acceptable (Sander, 2002), (iii) the formulation of hydraulic processes is very simple and should ideally be better constrained and (iv) organic matter only interacts with the soil through the action of increasing acidity so again missing out associated cation exchange and structural properties.However, many of these missing processes can be included in the model framework with relative ease once the method and relevant parameters are derived.For example, to best predict secondary mineral mineralogy the model could in future utilise a more complex chemical reaction module, for example by coupling with the PHREEQC geochemical programme (Parkhurst and Appelo, 1999).
Perhaps for this current study the largest limitation is that vegetation is prescribed in the model and productivity does not evolve with pedogenesis or nutrient availability.For this first attempt at pedogenesis modelling we prefer to keep the vegetation simple so that we can clearly identify how the vegetation may influence soil processes and development.
Introducing vegetation that responds to nutrient availability will require a more sophisticated vegetation module that is beyond the scope of this pedogenesis study.Instead, in the future this soil model may be coupled to an existing dynamic vegetation model.

Conclusions
This study has demonstrated that the soil evolution model presented is capable of reproducing realistic soil properties such as relative elemental losses, weathering depths, pH profiles, organic carbon content and soil-pore CO 2 concentrations.The model requires 20 parameters, of which at least 13 are easily assigned from literature, plus regional climate, bedrock data and simple thermodynamic constants to simulate soil genesis on a chosen parent material.The limited number of processes and the ease at which they can be both included and excluded from simulations makes the model behaviour easy to understand.This study has detailed how each of these model processes interacts with and influences the soil properties.
Comparisons of the model predictions with a Hawaiian soil chronosequences has highlighted the importance of vegetation in shaping soil profile evolution by increasing soil acidity and cycling nutrients.The good model agreement with the observations of Na, Mg and Ca which are less strongly cycled by vegetation, suggests that the model is realistically reproducing the other processes unrelated to nutrient cycling.These results lend confidence to the model's ability to quantify processes and feedbacks occurring during pedogenesis and to the valuable role it can play in understanding long-term biogeochemical cycles.assumed that solutions are ideal; thus mole fraction is equal to activity.This activity value is used to calculate the new ion concentration at thermodynamic equilibrium using the above procedure.
The model assumes that the behaviour of the elemental oxides depends only on their relative composition in the bedrock.However, these oxides are not usually present on their own, but are instead constituents of more complex silicate minerals.This will alter the solubility of the individual oxides and to account for this Kirkby (1977) proposed a correction term for the Gibbs free energy change of formation ( G f ) of each oxide.This correction term is determined by calculating the difference between the Gibbs free energy change of formation of the silicate minerals and the sum of the free energies of their constituent oxides.This difference is the formational free energy for the compound and is shared between the oxides to give the effective Gibbs free energy change of formation ( G f ).In this study a set of likely minerals is calculated from the weight percent of oxides in the parent rock and these are then used to find the correction factor.In order to determine the mineral assemblage of a rock from bulk chemical analysis a mineral norm is calculated.The norm is a set of idealised minerals that are calculated from the known composition of oxides in a rock.The method of calculating the minerals likely present is detailed below.

A2 Mineralogy and Gibbs correction factors
To calculate the likely minerals present in the parent material from the bulk chemical analyses of the rock, the CIPW (Cross, Idings, Pirsson and Washington) norm scheme is used.The CIPW method follows that of Hughes (1982).This method proceeds by expressing the oxides as molecular amounts and allocating the oxides to minerals in a step by step procedure.For example, all P 2 O 5 is used to make the mineral apatite which requires three times the amount of CaO, all TiO 2 is used to produce ilmenite using an equiv- alent amount of FeO, next all of K 2 O is used to make orthoclase and all Na 2 O is used to make albite unless there is not enough Al 2 O 3 , in which case the excess Na 2 O is used to make acmite.The method continues in this nature until all oxides are allocated to minerals.Some of the allocations are then revised depending upon the saturation or undersaturation of silica.FeO and Fe 2 O 3 are summed to give a total value for Fe and then a ratio of 0.1 Fe 3+ /total Fe is applied.This is because anomalously high Fe 2 O 3 contents can be recorded if the rock has undergone post-crystallisation oxidation (Hughes, 1982, pg. 97).The normative mineral assemblages obtained from the oxide compositions of Table A1 are shown in Table A2.This mineralogic configuration implies that the basalt is an alkali basalt where silica is undersaturated and nepheline is present.The correction factor for the G f for each oxide, λ, is found by comparing the free energies of these minerals with those of the constituent oxides and using the proportions of these minerals present in the parent material to find the mean λ for each oxide.For the model simulations of this study TiO 2 is assumed to be insoluble.

A3 Hargreaves equation
PET = 0.0023 × R a × (T mean + 17.8) × TD 0.5 , (A13) where PET is in units of mm month −1 .T mean is the monthly mean temperature ( • C), TD is the difference between the average monthly maximum and minimum temperatures ( is calculated for the 15th day of the month.The monthly value is calculated by multiplying this daily value by the number of days in the month.Equation (A14) estimates the extraterrestrial radiation using only latitude (φ) and the julian day (J ) (Kouwen, 2010).The proportion, p i , of oxide i remaining is then calculated from the original bedrock density and the loss of mass from leaching where m i (t = 0) is the mass of element i in the original parent material (ρ bedrock × wt % of i).

A5 Ionic diffusion
where D I is the diffusion coefficient of the ions which for current purposes we keep fixed for all elements.

A6 Bioturbation
where p i is the proportion of element i remaining in the profile at depth z and D (m 2 yr −1 ) is the diffusion coefficient.It is assumed that the mixing intensity will decline with depth due to the decrease in faunal activity with increasing soil depth (Humphreys and Field, 1998;Wilkinson and Humphreys, 2005;Johnson et al., 2014).In the model this takes the form of an exponential relationship: where D(0) is the diffusion coefficient at the soil surface and z b is the e-folding length scale for biological activity (m).
The boundary conditions at the top and bottom of the profile allow no mixing in or out so that ∂p i ∂z = 0, (A24) at z = 0 and z = zmax, where zmax is the total number of vertical layers in the model.

A7 Surface erosion
Erosion in the model acts by lowering the surface elevation.This lowering process shifts soil properties (or proportion of substance remaining, p) up the soil profile and thus

A8 Calculating chemical weathering intensity
To quantify element losses due to chemical weathering only, mass balance techniques can be applied to soil profiles (April et al., 1986;Riebe et al., 2004a, b).When soluble elements leave soil profiles the immobile elements become enriched compared to their concentrations in the parent material.Therefore measurements of immobile element enrichment in soil profiles can be exploited to reveal the extent of chemical weathering losses of other elements in the profile (e.g.Taylor and Blum, 1995).
To calculate the depletion or accumulation of an element relative to its concentration in the bedrock the soil/rock ratios of the element are normalised with those of a known inert element such as Zirconium (Zr) or Titanium (Ti) (Brimhall and Dietrich, 1987).This ensures that the differences in element concentrations between bedrock and soil is due to chemical weathering only and not because of changes in the soil bulk density or due to losses of other elements.
The weathering intensity of elements in the Hawaiian soils are calculated by normalising them relative to Nb: where τ is the mass transfer coefficient of element i (Brimhall and Dietrich, 1987), "p" stands for the protolith or parent material, "w" is the weathered material or soil and C is the concentration of element i. τ i = 1 indicates that the element is enriched at the same ratio as Nb and is therefore immobile, τ i = 0 indicates complete depletion of element i and τ i >1 represents relative enrichment of element i.

Figure 1 .
Figure 1.Diagram of the major processes, inputs and outputs of the soil profile model.

Figure 2 .
Figure 2. Vertical distributions of pH and the relative depletion/enrichment of the elemental model oxides after 20 kyr of soil development for five different model runs of increasing complexity.Values < 1 indicate a loss relative to the parent material, and values > 1 indicate relative accumulation.

Figure 3 .
Figure 3. Proportional mass loss of each model element relative to the amount in the parent material and mean pH for the top 50 cm of the soil profile after 20 kyrs of soil development for each of the model simulations (WL = weathering and leaching, E = erosion, B = bioturbation, C = organic carbon and N = nutrient cycling).

Figure 4 .
Figure 4. Proportion of original parent material remaining after 10 and 20 kyrs of model simulation for 3 different rates of surface erosion (T ).

Figure 7 .Figure 8 .
Figure 7. Observed vs. simulated oxide losses/gains across a mean annual precipitation gradient (P A ) on the 10ka Kona lava flow.Values < 1 indicate a loss relative to the parent material, and values > 1 indicate relative accumulation.

Figure 9 .
Figure 9. Observed vs. simulated oxide losses/gains across a mean annual precipitation gradient (P A ) on the 170 ka Hawi lava flow.Values < 1 indicate a loss relative to the parent material, and values > 1 indicate relative accumulation.

Figure 10 .
Figure 10.Observed vs. simulated oxide losses/gains across a mean annual precipitation gradient (P A ) on the 350 ka Pololu lava flow.Values < 1 indicate a loss relative to the parent material, and values > 1 indicate relative accumulation.The grey box shows the location of the pahoehoe flow.

Figure 11 .
Figure 11.Observed vs. simulated pH for 3 sites on each of the lava flows.

R
a = 15.392× d r (w s × sin φ × sin δ (A14)+ cos φ × cos δ × sin w s ) ,where d r is the relative distance between the earth and the sun, given byd r = 1 + 0.033 × cos 2π J 365 ,(A15)δ is the solar declination (radians) defined byδ = 0.4093 × sin 2π J 365 − 1.405 , (A16)and w s is the sunset hour angle (radians) given byw s = arccos(− tan φ × tan δ).(A17)A4 LeachingSee Fig.A1for a schematic of the leaching formulation.The amount of solute carried into a soil volume at depth z by percolating water is F (z)c i , where c i (g m −3 ) is the concentration of ion i in solution.The amount of solute lost from the volume element by diversion due to sub-surface flow is δF (z)c i and due to percolation outflow (F (z) + δF (z)) (c i + δc i ), thus−δ m i δz = [−(δF )c i − F c i + (F + δF )(c i + δc i )]δt, (A18)where m i is the mass change of oxide i at depth z during time δ t .Neglecting second order terms this reduces to p(z − δz, t + δt) = p(z, t),

Table A1 .
Oxide composition of a basalt and G f values for the rock oxides.