Modelling soil genesis in calcareous loess

Available online at www.sciencedirect.com

Geoderma 145 (2008) 462 – 479 www.elsevier.com/locate/geoderma

Modelling soil genesis in calcareous loess Peter A. Finke a,⁎, John L. Hutson b a

b

Department of Geology and Soil Science, University of Ghent, Krijgslaan 281, 9000 Ghent, Belgium School of Chemistry, Physics and Earth Sciences, Flinders University, GPO Box 2100, Adelaide 5001, Australia Received 5 February 2007; received in revised form 21 September 2007; accepted 29 January 2008 Available online 4 March 2008

Abstract The SoilGen1 model was developed to simulate soil development in calcareous loess at Holocene (15,000 BP–present) temporal extent. We used the LEACHC model as a core and added process formulations to describe the effect of vegetation and soil macrofauna on various soil properties. A limited calibration was done on the calcite solubility constant by comparing decarbonization rates in various leaching climates with values predicted by the metamodel of Egli and Fitze [Egli, M. and Fitze, P. 2001. Quantitative aspects of carbonate leaching of soils with differing ages and climates. Catena 46:35–62.]. Soil profiles of C and pH after 15,000 years compared well with measurements in a Belgian loess profile that was never under agriculture, taking the composition of the C-horizon as uniform parent material. Scenarios with and without agriculture, with varying degrees of bioturbation and for documented Holocene climatic and vegetation evolutions in Belgium and Hungary were simulated and compared. Results show a clear effect of bioturbation on soil development, especially in the continental climate evolution of Hungary. Indicators for the possible occurrence of clay migration and disturbance of lateglacial morphology were calculated and these also show the importance of bioturbation. © 2008 Elsevier B.V. All rights reserved. Keywords: Modelling; Soil genesis; Loess; Mechanistic model; Holocene

1. Introduction Many soils in the western European loess belt have in common (i) a parent material deposited in the late Weichsel period (up till about 15,000 BP), (ii) a mineralogical composition with some 10–15% calcium carbonate, some 12–18% clay minerals dominated by montmorillonite and most of the silt-sized fragments dominated by quartz minerals. As far as soil genesis in a climate with a precipitation surplus is concerned, there are often indications of decarbonization followed by clay migration. Clay migration may no longer be an active process in semi-natural environments, when pH drops to levels at which aluminium and iron enter the soil solution and keep the clays flocculated, while in agricultural environments liming may have favoured a restart of clay migration, often in combination with organic materials. In situations with a precipitation deficit, calcic horizons may be formed. Knowledge as summarized above, obtained by with

⁎ Corresponding author. Fax: +32 9 2644997. E-mail address: [email protected] (P.A. Finke). 0016-7061/$ - see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.geoderma.2008.01.017

mineralogical and soil chemical analysis and supported by soil morphological observations can be used to reconstruct soil genesis for particular pedons or toposequences (e.g., Brahy et al., 2000). Most of the underlying processes are now understood quite well and have been modelled, though not often at the temporal extent of soil formation. More generally, the paradigm of the so-called soil-forming factors that explain the soil that we observe is well established and accepted (Jenny, 1941; Dokuchaev, 1883), and most of the soilforming factors are used in mechanistic as well as empiric soil models to satisfy initial or boundary conditions. Nevertheless, only few attempts to model soil formation at large temporal extents have been reported (Schaetzl and Schwenner, 2006), and mostly predict sets of soil variables relevant to subsets of processes in the soil like soil acidification (e.g., Kros et al., 1999) and C-sequestration (e.g., Parton et al., 1993). Soil formation models, especially mechanistic models, could be useful temporal interpolators between parent material and present-day soil (e.g. to assess the agricultural carrying capacity in an archaeological context, the effect of the introduction of agriculture on soil development), temporal extrapolators (e.g. in the context of

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

climate change), and could provide benchmark points for soilscape reconstruction and digital soil mapping (McBratney et al., 2003). Mechanistic models for soil genesis are not easily operationalized because of the following reasons: (i) The temporal variation of boundary conditions such as precipitation and air temperature is not easy to reconstruct. Recently however, this situation has – in Europe – much improved now that local and regional climate reconstructions have become available (Petit et al., 1999; Seppä and Birks, 2001; Davis et al., 2003). (ii) With increasing temporal extent an increasing number of processes must be described in the model to evaluate the effects of (e.g.) lateglacial temperatures in soils of the temperate regions; of vegetation change and associated change in organic C; of changes in the vertical distribution of solid phase components. This increases the models' complexity and associated effort to estimate parameters and may decrease the models' stability. (iii) Limited experience exists with the mechanistic modelling of some processes relevant for soil formation such as bioturbation and clay migration. (iv) Calibration of mechanistic models of soil formation is difficult because often only the current stage of soil development can be measured quantitatively. In some cases, this situation is more favourable when chronosequences in the same parent material are found. Nevertheless, model calibration using state variables may be less feasible than calibrating selected process rate variables or some temporal integral of these. This situation is less critical when the model has a positive verification status in similar soil environments. (v) Computational demands: the runtime of such model may be high. This paper reports on modelling soil formation in a calcareous loess parent material by confronting the above problems. The

463

research objectives were: (i) to develop a process-based model based on present-day knowledge that is able to reconstruct (postglacial) soil development in calcareous loess; (ii) to calibrate, test and apply the model in two locations in the loess belt of western and central Europe with contrasting climatic evolutions from 15,000 BP to present. The first objective leads to the use of a model of documented high quality for water and solute transport as the core for the soil genesis model. The LEACHM model (Hutson and Wagenet, 1992), more specifically the model variant LEACHC for the transport of inorganic ions in soils, was chosen because of its documented quality, wide application range and numerical stability. A substantial number of functions must however still be added to this core to cover the wider range of pedogenetical processes. One of these is the cycling of C, for which the concepts of the RothC 26.3 model (Jenkinson and Coleman, 1994) were implemented. The resulting model, which is an integration of LEACHC translated from Fortran to the Delphi™ Pascal platform and newly developed codes, will be referred to henceforth as SoilGen1. 2. The SoilGen1 model 2.1. Model structure Table 1 summarizes how the Jenny (1941) factors of soil formation are linked to the model and what existing model (LEACHC), model description (RothC 26.3) or new functionality was implemented to model the soil processes influenced by the Jenny factors. Reference is made to Hutson (2003a,b) and to Coleman and Jenkinson (2005) for detailed descriptions of the LEACHC and RothC 26.3 models. Although the temporal extent is large (15,000 years), the timesteps within the model may be small, depending on the dynamics of the processes. Most processes are modelled on the subday timescale. Fig. 1 gives a process order diagram within one year. The flow of water and solutes is computed with a maximum timestep of 0.05 day, and this is reduced in case of large fluxes. CO2-diffusion is computed with timesteps of 0.1 day, which is

Table 1 Implementation of soil-forming factors in SoilGen1 Environmental factor Climate

Organisms

Relief Parent material

Time

Temperature Precipitation: water Precipitation: solutes Evaporation Vegetation

Fauna Human influence Slope Variants of T, P, E Texture Mineralogy CaCO3 + CaSO4 + salts Change of all BCs

Link to model

Modelled processes

Source of code

BC BC BC BC BC SIM SIM BC BC BC IC BCs IC + SIM IC IC + SIM BC

Heat flow Water flow Solute flow Evapotranspiration C-cycle CO2-production and diffusion Selective cation uptake/release Root distribution Bioturbation Fertilization Runoff Heat/water/solute flow Effects of dissolution/ precipitation, bioturbation and C-cycling Cation exchange Chemical equilibria Annual update of all BCs

LEACHM LEACHC LEACHC LEACHC/modified After RothC 26.3 New New New New LEACHC LEACHC LEACHC New LEACHC LEACHC New

BC, IC = boundary, initial condition, SIM = simulated.

464

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

reduced in case of numerical instability. The recalculation of chemical equilibria is done once daily or every 80 water flow timesteps, whichever value is smaller. The physical redistribution of soil matter due to dissolution/precipitation and bioturbation is calculated once per year since these changes are slow. 2.2. Model components 2.2.1. Flow of water, solutes, heat and CO2 The flows of water, solutes and heat are simulated using the LEACHC code. Below follows a brief characterization of the model. For details reference is made to Hutson (2003a). Water flow is modelled by Richards' equation for transient vertical flow:
 Ah A AH C ð hÞ ¼ K ð hÞ  U ð z; t Þ ð1Þ At Az Az where C(θ) is the differential water capacity ∂θ / ∂h, θ is the volumetric water content (m3/m3), h is soil water pressure head (Pa 10), K(θ) is hydraulic conductivity (m 10− 3/d), H is hydraulic head (Pa 10) and U(z,t) is a sink term representing water lost at depth z and time t by transpiration. LEACHC solves the Richards' equation by application of finite differencing techniques. Steps taken are the rearranging and subsequent application of the Crank–Nicholson implicit method to estimate pressure heads at the end of the timestep and by application of a Gaussian elimination method (Thomas algorithm) to simultaneously solve the equations at all nodes in the soil profile. Note that all node distances were set to 5 cm and that total profile depth was set to 150 cm for all simulations. The upper and lower nodes are outside the soil and are used to impose the upper boundary conditions for evaporation and infiltration and the lower boundary condition. The lower boundary condition is chosen by the model to cope with strong variations in the precipitation and evaporation regime over the last 15,000 years: zero flux in case of an annual precipitation deficit and free drainage in case of an annual precipitation surplus. The modified Campbell equations (Campbell et al., 1977; Hutson, 2003a,b), are used to obtain values for the differential water capacity and hydraulic conductivity from values of θ or h. To account for the effect of temporarily frozen soils, the hydraulic conductivity was recalculated for soil temperatures below 0 °C using an impedance factor Ω (Lundin, 1989): K = 10−Ω K. Ω was assigned the value 4. Heat flow and temperature distribution are modelled by (Hutson, 2003a; Tillotson et al., 1980):   AT A Kt ðhÞ AT ¼ d ð2Þ At Az b Az where T is temperature (°C), Kt(θ) is thermal conductivity (J m− 1 s− 1 °C− 1) calculated at θ using the method presented by Wierenga et al. (1969) and β is the volumetric heat capacity determined from β = ρsCs + θCwρw with ρs and ρw the bulk densities of solids and water (1000 kg m− 3) respectively, Cs the gravimetric heat capacity of solids (840 J kg− 1 °C− 1) and Cw the gravimetric heat capacity of water (4200 J kg− 1 °C− 1). Eq. (2) is solved for all profile nodes using an implicit centraldifference scheme with a Gaussian elimination method. The upper boundary condition is satisfied by a sinusoidal daily air

Fig. 1. Process flow (arrows) and temporal scales (solid boxes or dots) of the subprocesses in SoilGen1 in each 1 year. Open boxes indicate groups of processes that are repeated at daily of annual extent.

temperature fluctuation derived from (input) daily averages and amplitudes of temperature. The lower boundary condition is a heat reservoir: the deepest compartment is assigned a thickness of 2 m and there is a zero heat flux bottom boundary condition. The flow of soluble matter is simulated using a finite difference approximation to the convection–dispersion equation (CDE) for each soluble compound: AðhC Þ A ¼ ½hDðh;qÞ  qC FU At Az

ð3Þ

where C is the solute concentration (kg m− 3), D(θ,q) is the dispersion coefficient (mm2 d− 1) being the combined effect of mechanical dispersion and aqueous diffusion, q is the water flux (mm d− 1) and Φ is a source or sink term (kg m− 3 d− 1) representing the plant uptake or release by mineralization of organic matter. Note that solute concentrations depend on chemical equilibria and

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

partitioning between the exchange and solute phase (cf). A centraldifference Cranck–Nicholson approach is applied while solving Eq. (3) to avoid numerical dispersion. Upper boundary conditions allow for surface infiltration, evaporation and zero flux, while the lower boundary condition depends on that for water flow: zero concentration for zero flux and constant concentration for free drainage. The flow of CO2 is assumed to be diffusive and is simulated by an explicit numerical solution to the gas regime equation: ed

Ac A2 ¼ DðT Þgs d 2 þ Pð z; t Þ At Az

ð4Þ

where ε is the air-filled porosity, c is the CO2-concentration (partial pressure) in the soil air, P(z,t) is the CO2 production in each soil compartment and D(T)gs is the gas diffusion coefficient in soil (m2 s− 1), estimated by (Moldrup et al., 2000):  2þ3=b  3  e DðT Þgs ¼ DðT Þ0 d 2e100 þ 0:04e100 d ð5Þ e100 with ε100 is the air-filled porosity at −100 cm pressure head, b is the Campbell soil water retention parameter also used for water flow calculations and D(T)0 is the gas diffusion coefficient in free air obtained by (assuming a constant pressure of 101.3 kPa):   T þ 273:16 1:75 ð6Þ DðT Þ0 ¼ 1:39  105 d 273:16 with T in °C. 2.2.2. Plant related processes SoilGen1 distinguishes 4 vegetation types (grass/scrub, conifers, deciduous wood and agriculture), each of which is

465

characterized by: (i) a rooting density function, (ii) preferential cation uptake controlled by fixed cation ratios in the plant, (iii) C-decomposition rates, (iv) the fraction of dead C that enters the soil system as leaf and root litter, (v) the distribution of litter input over the year and (vi) a partitioning coefficient to fraction the fresh litter in resistant and decomposable components. Furthermore, the yearly amount of produced C-litter is input to the model (and normally depends on the vegetation). The belowground C-cycle (visualized in Fig. 2) is modelled according to the concepts of the RothC 26.3 model (Jenkinson and Coleman, 1994). Dead plant material, split in an ectorganic (leaf litter) and endorganic (root litter) part (Kononova, 1975), is divided into a resistant and a decomposable fraction according to a vegetation-dependent ratio. Both fractions decompose into biomass, humus and CO2 at rates that are determined by the fraction that is decomposing, soil temperature, soil moisture deficit, soil cover fraction and the time increment. Biomass and humus continue to decompose into biomass, humus and CO2 in next timesteps. For a detailed description and values of rate factors reference is made to Coleman and Jenkinson (2005). The C-submodel is applied at daily intervals, where the produced CO2 enters the gas regime equation (Eq. (4)). The CO2 profile at the end of each day gives the pCO2 values for the chemical equilibria for this day. Uptake of the cations Ca, Mg, K and Na by vegetation occurs via the transpiration stream and is forced to reflect the relative proportions of those elements measured in the plant (stable vegetation assumption). Each vegetation is characterized by a content of the basic cations Ca, Mg, K and Na. Typical contents for agriculture (Barley) were obtained from Wyszkowski et al. (2006), for grass/scrub from Thompson et al. (1997), and for coniferous and deciduous wood from Navrátil (2003). These

Fig. 2. C-cycling in SoilGen1 (based on Coleman and Jenkinson, 2005). Greyshaded boxes indicate pools, rounded boxes indicate processes and the rectangular white box is for conceptualization only.

466

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

Table 2 Relative concentrations of Ca, Mg, K and Na in 4 vegetation types Vegetation

Coniferous forest Deciduous forest Grass/scrub Agriculture

Source

Navrátil (2003) Navrátil (2003) Thompson et al. (1997) Wyszkowski et al. (2006)

Relative concentration (−) Ca

Mg

Na

K

0.770 0.483 0.271 0.134

0.131 0.227 0.165 0.076

0.004 0.025 0.004 0.013

0.095 0.266 0.560 0.776

contents are rescaled to 4 cation mass fractions summing up to 1 (Table 2). The calculation of cation and anion uptake per soil compartment then proceeds as follows: 1. A general cation uptake factor is calculated as the minimal occurring ratio between one of the 4 cations actually present in the solution in the soil compartment and the associated cation mass fraction. 2. An uptake fraction for each cation in the soil compartment is calculated by the multiplication of the general cation uptake factor and the cation-specific mass fraction. This fraction is applied to this cation in solution in the soil compartment to calculate actual plant uptake. If one of the cations has 0 concentration, non-preferential uptake is assumed. 3. The uptake charge from the soil compartment in the transpiration stream is calculated, and it is checked if the sum of the anion charges of Cl, SO4, CO3, HCO3 in solution in the soil compartment can balance this charge. There is no preferential uptake of any anion. If there is no charge balance, a correction ratio to the cation uptake is calculated satisfying 0-charge of total uptake.

The absolute uptake in one timestep is thus limited by the element concentrations in the rooted soil compartments, by the transpiration flux and by a charge balance condition. Cations are stored in the same plant biomass pools distinguished in the C-cycle, and those cations entering the mineralized (CO2) pool actually enter the soil solution either at the soil surface (leaf litter) or in the rooted compartments (root turnover). 2.2.3. Chemical processes The solution phase is brought in equilibrium with the precipitated and exchange phases by satisfying the following thermodynamic equilibria: (i) Henry's Law constant for CO2, (ii) the dissociation constant of H2CO3, (iii) the dissociation constant of water, (iv) the solubility constants of gypsum and calcite, (v) ion pair stability constants for the species named in Fig. 3, (vi) Gapon selectivity constants for the exchange/ solution phase equilibria for Ca–Mg–Na–K. Calculation of the equilibrium distribution is done via an iterative procedure: (i) for the initial values of the pH and elements in solution, the ionic strength is calculated, which is used to calculate activity coefficients for all ion species and equilibrium constants using the Davies' relationship (Stumm and Morgan, 1970); (ii) the levels of HCO3−, CO32−, CaSO4 and CaCO3 are calculated using pH and pCO2 (which is daily adjusted); (iii) the cation levels are partitioned between the solution and the exchange phases according to the Gapon equations; (iv) steps (i) to (iii) are repeated until the result is stable;

Fig. 3. Inorganic chemical processes (arrows), boundary conditions, sources and sinks (white boxes) and chemical phases (greyshaded) in SoilGen1 (adapted after Hutson, 2003).

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

(v) with a stable result, the charge balance is checked, and if an unbalance exists, pH is adjusted and the procedure restarts at step (i). For a detailed description of the calculations involved, reference is made to Hutson (2003a). 2.2.4. Processes causing redistribution of soil phases The processes considered in SoilGen1 able to change the distribution of solid phase and liquid phase components are bioturbation, tillage and mass changes due to accumulation or mineralization of OC and dissolution or precipitation of calcite and gypsum. Erosion and sedimentation are not considered. Central assumption in SoilGen1 is that the volume of all soil compartments is constant over time, although in reality a volume of soil may lose porosity (collapse) due to the removal of calcite, or may gain porosity due to biological activity. Thus, in the terminology of Brimhall and Dietrich (1987) a strain equal to 0 was assumed. The errors thus introduced do not affect calculated mass percentages, but may have some influence on soil physical characteristics. Redistribution of mass by bioturbation is done in 2 steps: (i) The percentage of the mass subject to vertical redistribution by soil meso- and macrofauna in each compartment is determined. Currently, this percentage is input and was made to vary over the Holocene time extent with the vegetation, climate and soil depth. Whole-soil values of bioturbation for different vegetations from Gobat et al. (1998, p.122) were taken as reference. Bioturbation mass percentages of all mineral pools (clay, silt, sand, including calcite and gypsum), soil water and elements in solution and OC-pools are put in the vertical mixing pool. In this pool, values are averaged and masses are reassigned to the bioturbated soil compartments according to the bioturbation mass percentages. (ii) For each bioturbated compartment, the bioturbated mass is mixed with the non-bioturbated mass to obtain one set of soil properties per soil compartment. An additional effect of bioturbation is the increasing gas diffusion. Singer et al. (2001) measured a CO2-diffusion of 5 ×10− 4 cm2/s in the absence and 4.45 × 10− 3 cm2/s in the presence of earthworms. This factor 8.9 increase of D(T)gs (Eq. (5)) is applied to those compartments with bioturbation. Mass changes, as induced by bioturbation but also by changes in OC, calcite or gypsum, in all compartments are used to (annually) recalculate mass percentages of the solid phase components OC, Clay, Silt and Sand. Subsequently, bulk density and porosity are recalculated and pedotransfer functions (Rawls and Brakensiek (1985) are applied to recalculate the Campbell parameters (Campbell et al., 1977). 2.3. Model calibration and tests The usual way to calibrate a model is to estimate to what model constants the output of selected variables is most sensitive, and then to tune these constants so that measurements

467

Table 3 Initial conditions and some constants for all soil layers in all simulated scenarios at 15,000 BP Type

Variable/constant

Value

Physical

Clay (a) Silt (a) OC (a) Temperature PD clay PD silt + sand PD–OM Bulk density Water content Ca2+ solution (a,d) Mg2+ solution (a,d) Na+ solution (a,d) K+ solution (a,d) Cl− solution (a,d) SO2− 4 solution (a,d) Alkalinity (a) Ca-exchange (a,d) Mg-exchange (a,d) Na-exchange (a,d) K-exchange (a,d) CEC (a) pCO2 CaCO3 (a) CaSO4 (a) Wilting point Min. root water potential Max ratio act/pot T Root resistance Plant cover Slope log(KSO) CaCO3 (e) log(KSO) CaSO4 KG Mg/Ca (b) KG Na/Ca(b) KG K/Ca (b) KG Na/Mg (b) KG K/Mg (b) KG K/Na (b) Ca–rainwater (c) Mg–rainwater (c) Na–rainwater (c) K–rainwater (c) Cl–rainwater (c) SO4–rainwater (c) Alkalinity–rainwater (c)

12.6% 73.4% 0.12% −5 °C 2.65 kg/dm3 2.65 kg/dm3 1.10 kg/dm3 1.30 kg/dm3 0.30 dm3/dm3 15.91 mmol/dm3 1.17 mmol/dm3 0.07 mmol/dm3 0.07 mmol/dm3 34.0 mmol/dm3 0 mmol/dm3 0.30 mmol/dm3 68.59 mmol+/kg 13.01 mmol+/kg 0.70 mmol+/kg 1.50 mmol+/kg 83.80 mmol+/kg 0.1 kPa 12% 0% −1500 kPa −3000 kPa 1.1 1.05 100% 0% −8.46 −4.61 0.6982 14.9 30.2 21.37 43.27 2.02 0.005 mmol/dm3 0.010 mmol/dm3 0.140 mmol/dm3 0.010 mmol/dm3 0.150 mmol/dm3 0.035 mmol/dm3 −0.040 mmol/dm3

Chemical

Plant

Land Chemical

PD = particle density, T = transpiration, KG is Gapon exchange constant, KSO is solubility constant. Sources: (a) Van Ranst, 1981, (b) De Vries and Posch, 2003, (c) Monitoring data 2005, (d) extractable cations, (e) calibration. Cations (a,d) were distributed over exchange and solution phase using (b). Non-referred values are own estimates.

are most accurately reproduced. When simulating pedogenesis, the only measurements available are at the start (parent material) and at the end of the simulation period, unless sampling locations in the same parent material can be found where soil formation has started later (chronosequences). The loess deposit in which soil formation is simulated is the most recent loess deposit in Europe, and restarts of soil formation after soil erosion cover only very brief periods since soil erosion in loess is most often associated with agriculture and thus a more or less continuous process. This implies that any calibration aiming at a

468

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

Fig. 4. January and July temperatures for Belgium and Hungary (after Davis et al., 2003).

comparison between measured and simulated finite state will be based on few data and will also be time-consuming because of the runtime of the model. We therefore decided: (i) To rely on the verification status of the model components LEACHC and RothC on short timescales, as documented in (a.o.) Smith et al. (1997), Addiscott and Wagenet (1985), Jalali and Rowell (2003), Dann et al. (2006) and Jabro et al. (2006); (ii) To use finite state measurements for verification only; (iii) To calibrate the model if possible on process rates rather than system states; (iv) To obtain an impression of the sensitivity of the model to added processes like bioturbation by defining simulation scenarios and comparing the results.

One crucial process in the model is the rate of decarbonization, since it will determine the development of the cation content in the exchange and solution phases, pH, porosity and physical characteristics. Egli and Fitze (2001) analysed decarbonization rates for various soils in Europe and developed a metamodel to predict decarbonization rates (mol m− 2 y− 1) as a function of yearly soil water flux and soil texture. The decarbonization predicted by such metamodel can be compared to decarbonization predicted by SoilGen1 if the annual leaching flux, calcite content and soil texture are known. As the decarbonization speed in SoilGen1 is largely determined by the solubility constant (KSO) of CaCO3, we calibrated this constant. We used the silt loam texture and calcite content of the parent material (Table 3) at a high annual percolation for the Belgian climate in the Holocene (472 mm) to derive the number of years needed to decarbonize

Fig. 5. Annual precipitation (after Davis, unpublished data) and potential evapotranspiration (Hargreaves method) for Belgium and Hungary.

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

469

Fig. 6. Dominant vegetation (simplified after Verbruggen et al., 1996) and assumed bioturbation for Belgium. DE = deciduous wood, G/S = grass/scrubland, CO = coniferous wood.

1100 mm of parent material with the Egli and Fitze (2001) metamodel, and compared this with SoilGen1 outputs for various KSO. We assumed an annual bioturbation of 30 ton ha− 1 and a deciduous forest for vegetation. The value of KSO giving the best match between the decarbonization times of both models was used henceforth. Thereafter, SoilGen1 with the chosen value for KSO was run on the same input data, but for different values of the annual percolation (considered representative (Fig. 5) for the Holocene climate variation in Belgium). To check model performance with the chosen KSO, the decarbonization time of 1100 mm of parent material was compared to the result obtained with the Egli and Fitze (2001) metamodel.

3. Exogenous model inputs on the Holocene timescale To evaluate the effect of the soil-forming factors on soil formation with the model, time series of these factors must be generated as input. We assumed constant relief (plateau position) and parent material (no erosion, no sedimentation) over the Holocene, although the model is capable to simulate soil formation in layered parent materials and can also handle erosion and sedimentation to some extent by addition or removal of compartments at the top of the simulated profile. Time series over 15,000 years had to be estimated for the climate variable temperature, precipitation and potential evapotranspiration and

Fig. 7. Dominant vegetation (simplified after Rybnícknova and Rybnicek, 1996) and assumed bioturbation scenarios for the Hungarian plain. G/S = grass/scrubland, CO = coniferous wood.

470

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

Table 4 Characterization of simulated soil formation scenarios Climate

Natural development

Agriculture

Bioturbation level (ton soil ha− 1 y− 1)

Bioturbation level

None (0)

Normal (8–16)

Belgium BEL Nat NoBio Hungary HUN Nat NoBio

Normal (8–16)

High (16–32)

BEL Nat Bio BEL Agr Bio HUN Nat Bio HUN Nat HiBio HUN Agr Bio

for the variables describing the impacts of soil organisms (vegetation, bioturbation) and men (absence/presence of agriculture and fertilization, more specifically, liming). These time series served as boundary conditions to the model. 3.1. Climate Temperature data are used by the heat flow subroutine to calculate soil temperatures. Soil temperatures directly influence decomposition rates of organic matter and CO2 production rates and thus indirectly the decarbonization. Freezing temperatures directly affect the hydraulic conductivity of the soil and thus the moisture distribution in the profile and leaching and possibly also surface runoff. Davis et al. (2003) published time series of area-averaged air temperature anomalies for January and July for 100-year intervals of the most recent 12,000 years, based on quantitative analysis of pollen data. The temperature anomaly series for central-western Europe and central-eastern Europe (Fig. 4) were used in this study. Temperature anomaly data from the Vostok Ice Core Data set (Petit et al., 1999) for 15,000 BP– 12,000 BP were combined with the Davis et al. (2003) temperature anomaly data by application of a timeshift and temperature anomaly shift so that the temperature anomalies during and after the Alleröd and Younger Dryas maximally coincided. This resulted in an extension of the year-average temperature anomalies for 15,000–12,000 BP. The thus obtained temperature anomaly time series 15,000BP–present for 2 regions and for July and January were converted into actual temperature series of July and January temperatures by

adding current data from Uccle (Belgium) and Gödöllő (Hungary). The resulting values for each 100 years (Fig. 4) were linearly interpolated to obtain yearly values, and subsequently used to scale a standard time series of weekly temperature data so that January and July values were satisfied. Scaling was done by calculating the difference of the standard time series and the climatic time series for the months January and July (annually), interpolating this difference for the other months and adding the monthly differences to the standard series. Within-daily temperature fluctuations are simulated by a sinusoidal relationship with time, assuming a daily amplitude from this standard time series. Precipitation data and potential evapotranspiration data are used to calculate water flow and evapotranspiration, and thus directly influence the leaching and uptake of solutes. There are many indirect influences such as the CO2 production rate and the partial CO2-pressure in the soil. Annual precipitation data anomalies for both European regions for the most recent 12,000 years were supplied by Davis (unpublished data). These data were converted to actual precipitation values by adding current annual values for Uccle and Gödöllő (Fig. 5). The time series were extended to 15,000 BP by assuming minimal values for precipitation in the Lateglacial and Younger Dryas periods as reported by Liu et al. (1995) and guessing values for Alleröd, Middle Dryas and Bölling periods. Yearly values were obtained from the 100-year values by linear interpolation. These yearly values were used to scale a standard series of daily precipitation data by a multiplier. The method of Hargreaves and Samani (1985) was applied to calculate daily values for ET0 for 51° NL using the current monthly average temperatures and daily temperature ranges of Uccle. The yearly sum of ET0 was calculated and a correction factor was obtained by comparison with measurements from Uccle. Assuming that daily temperature ranges would not change with changes in annual average temperature, the calibrated model was applied to temperatures corresponding to annual averages between −6.5 and 15 °C. The obtained ET0 were used to obtain a regression model of the form ET0 = a ⁎ Tyr + b. This regression model was used to estimate ET0 using the earlier obtained temperature data (Fig. 5).

Fig. 8. Calibration (left) and verification (right) of the solubility constant for calcite at annual precipitation surplus of 472 mm (left) and 247, 307, 412, 472 and 516 mm (right) of SoilGen1 against a metamodel from Egli and Fitze (2001).

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479 Table 5 Output variables Type

Variable

Dimension

Chemical

Alkalinity Ca2+ solution Mg2+ solution Na+ solution K+ solution Cl− solution SO2− 4 solution CO2− 3 solution HCO−3 solution EC pH Ca-exchange Mg-exchange Na-exchange K-exchange ESP SAR pCO2 CaCO3 CaSO4 Clay Clay Sand Silt RHO Porosity Theta Potential ET Flux Temperature 31–12 Average temperature OC DPM b RPM b Biomass b Humus b IOM b Roots

mmol− dm− 3 mmol dm− 3 a mmol dm− 3 a mmol dm− 3 a mmol dm− 3 a mmol dm− 3 a mmol dm− 3 a mmol dm− 3 a mmol dm− 3 a mS m− 1

Physical

Plant related

mmol+ kg− 1 soil mmol+ kg− 1 soil mmol+ kg− 1 soil mmol+ kg− 1 soil % – bar mass fraction mass fraction mass% fine earth volume% mass% fine earth mass% fine earth kg dm− 3 volume fraction cm3 cm− 3 kPa mm mm °C °C mass% solid fraction ton C ha− 1 ton C ha− 1 ton C ha− 1 ton C ha− 1 ton C ha− 1 Fraction

All variables per soil compartment (5 cm, 0–150 cm) and at 31–12 of each year unless indicated otherwise. DPM = decomposable plant material, RPM = resistant plant material, IOM = inert organic matter, OC = organic carbon. a Also as mmol m− 2 in ectorganic matter. b Also in ectorganic matter.

3.2. Organisms Vegetation data are used in the model to derive the root density distribution and the preferential uptake of cations, to calculate transpiration and to calculate daily input of litter at the soil surface and in the soil. Indirectly, the organic pools determine the CO2-production and the release of ions via mineralization. Vegetation data for Belgium and Hungary were obtained by generalizing vegetations from pollen data (Verbruggen et al., 1996; Rybnícknova and Rybnicek, 1996) into 4 vegetation types: grass/scrubland, coniferous forest, deciduous forest and agriculture (Figs. 6 and 7). Bioturbation data are used to determine the degree of mixing of the topsoil compartments and to increase CO2-diffusion. Indicative values from Gobat et al. (1998, p.122) for different

471

vegetations were used to define the depth of maximal bioturbation, maximum depth of bioturbation and the mass fractions of the solid phase mixed at the soil surface, at depth of maximal bioturbation and at the maximum depth of bioturbation in such way, that the vegetation-dependent indicative bioturbation values in ton ha− 1 y− 1 were reproduced (Figs. 6 and 7). Influence of men input to the model as the temporal extent of the vegetation type “agriculture” and as liming. The starting date of agriculture was taken from the pollen diagrams (Verbruggen et al., 1996; Rybnícknova and Rybnicek, 1996) and set to 2560 BP for Belgium and 2960 BP for Hungary. Agriculture in both regions has a much older history, but was probably not continuous in time and over large areas before these years. According to Plinius the Young liming was done in Roman times once every 10 years by application of about 1 ton marl/ac. Assuming 75% CaCO3 in marl, this gives 1.875 mol CaCO3 m− 2 applied every 10 years. Current maintenance liming in Belgium is equivalent to 0.276 mol CaCO3 m− 2 y− 1. We assumed the 10-year cycle of liming from the start of agriculture to 60 years BP and maintenance liming afterwards. 3.3. Scenarios of soil formation Seven scenarios of 15,000 years of soil formation were simulated (Table 4). These scenarios varied in the evolution of the climate and vegetation, the intensity of bioturbation and the presence or absence of agriculture. All scenarios started with the same parent material: frozen calcareous loess containing a limited amount of (inert) organic matter (Table 3). 4. Results and discussion 4.1. Calibration and verification Fig. 8 shows that at a log(KSO) of −8.46, the time to decarbonize 1100 mm of loess of both the SoilGen1 model and the Egli and Fitze (2001) metamodel is equal to 1297 years. Application of this value to a range of annual precipitation surpluses between 247 and 516 mm yielded decarbonization times that compare well (RMSE of 67 years) with values obtained by using the metamodel, although at low precipitation surpluses SoilGen1 decarbonization values are systematically lower (Mean error of 45 years). The C-content and the pH of the soil after 15,000 years were compared with measured values from the Zonian forest, Belgium (Fig. 9). This location has never been under agriculture and thus the soils represent natural soil development. Both the soil pH as the carbon content compare well with model predictions. Measured organic matter content close to the soil surface is higher than simulated, and the simulated carbon content is more gradually decreasing with depth than the measured profile. Partly this can be attributed by the sampling depths and the smoothing effect of the 5 cm compartments in the Soilgen1 model, but possibly this also indicates an overestimation of (recent) bioturbation (input to the model). The simulated low pH values are near the limit of the application domain of the model, since aluminium chemistry

472

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

Fig. 9. Verification of predicted C% and pH after 15,000 years (undisturbed loess, Belgium).

and mineral weathering are not accounted for and thus the buffering mechanism in the pH-range near 4 is absent. 4.2. Scenario output and postprocessing For each scenario, the soil variables mentioned in Table 5 were calculated per year and per soil compartment. A selection of the output variables was depicted in time–depth diagrams (Figs. 10–16) and is discussed hereafter. The organic carbon content (Fig. 10) is related to the vegetation (model input) and clearly shows the effect of the introduction of agriculture (BEL Agr Bio versus BEL Nat Bio scenario, Table 4). Some variations occur over short periods and are related to discrete changes in vegetation. Highest values are found in the Bölling/Alleröd periods, when the production of C

strongly increased to present-day levels and decomposition was less than present because of low winter temperatures. In reality, soil OC levels may have been higher in these periods because of the stabilisation of OC by Ca that was not yet leached. This process is however not included in the model. The levels of organic carbon at the soils surface (ectorganic matter, Fig. 11) are highest during the Bölling/Alleröd periods as well. A strong decrease can be observed for the Younger Dryas period due to the herbaceous vegetation and in the Holocene the level of ectorganic matter shows an initial increase and subsequently stabilisation. The calcite content (Fig. 12) for the BEL scenario reflects the leaching climate: decarbonization of the upper 120 cm is complete after less than 2000 years. The degree of decarbonization in the HUN scenarios depends on the degree of bioturbation

Fig. 10. Simulated C-content.

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

and if liming is applied (HUN Agr Bio). All HUN scenarios show a loss of calcite in the topsoil and the formation of an horizon of CaCO3 enrichment at depth near 75 cm. These meet the WRB diagnostic criteria for a calcic horizon (IUSS Working Group WRB, 2006: more than 15% CaCO3, an increase of more than 5% CaCO3 relative to the underlying layer and a thickness 15 cm or more). The scenario without bioturbation leads to a more pronounced decarbonization and associated build-up of a calcic horizon than the HUN Nat Bio, HUN Nat HiBio and HUN Agr Bio scenarios. Especially in the period 10,000 BP– 6000 BP, bioturbation brings sufficient calcite to the surface to keep pH high and clearly slows down decarbonization for thousands of years. Bioturbation slows down decarbonization in BEL scenarios as well, but only for tens of years. Changes in the dry bulk density (Fig. 13) mainly reflect the changes in calcite content. The BEL Nat Bio scenario shows the effect of decarbonization on bulk density. The HUN Nat Bio and HUN Agr Bio scenarios show the effect of liming and decarbonization in the topsoil. The more pronounced calcic horizon that develops in the HUN Nat Bio scenario is visible as well. These statements are subject to the assumption that these soils developed in the absence of strain (Brimhall and Dietrich, 1987).

473

The effect of vegetation type on the distribution of exchangeable K (Fig. 14) is clearly seen. These cations concentrate in the topsoil as long as deciduous wood vegetation persists (HUN Nat Bio), but are quickly leached after agriculture starts (HUN Agr Bio). This pattern occurs as well for exchangeable Ca and Mg in the HUN and BEL scenarios (not shown), but not for Na. It should be noted that values of exchangeable cations Ca2+, Mg2+, Na+ and K+ are indicative and probably overestimated at low soil pH as exchange acidity (Al3+ and H+) is not estimated by the model. In addition to the mechanistic description of the processes identified in Table 1, two indicators for the occurrence of soil processes and soil properties were obtained via postprocessing: (i) an indicator variable for the possible occurrence of clay migration; (ii) an indicator for the preservation of frost-related soil morphology, as often visible in Pleniglacial deposits such as loess (Van Vliet-Lanoë et al., 1992). The state of progress in the deterministic description of clay migration does not yet allow (DeNovio et al., 2004) inclusion of this process in a mechanistic model. The possible occurrence of clay migration was therefore modelled by a clay migration indicator ICM, using an aggregate stability diagram (Van Breemen

Fig. 11. Simulated ectorganic C. L = decomposable organic matter, F = resistant organic matter, H = humus + biomass + inert organic matter.

474

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

Fig. 12. Simulated CaCO3 content.

and Buurman, 2002: p.195) which was uncertainty in the threshold values: 8 0 8 pH > > < pH  4 8 pH ICMpH ¼ 2  ð7:5  pHÞ 8 pH > > : 1 8 pH

fuzzified to account for g a a a

9 ½4;7:5 > > = ½4;5 ½7;7:5 > > ; ½5;7

ð7Þ

where pH is calculated for each soil compartment and year. A value of ICM of 1 indicates likely, a value of 0 indicates unlikely clay dispersion. It should be noted that a high likeliness of clay dispersion does not necessarily indicate massive clay transport. Other factors such as the presence of air–water interfaces and variations in pore water velocity are of major importance as well (DeNovio et al., 2004). Fig. 15 depicts the development of ICM in loess soils under natural and agricultural conditions for Belgian and Hungarian climate developments. The effect of soil liming in the agriculture scenarios is clearly visible: in HUN Agr Bio the

time-window of possible clay migration closes after liming takes place because pH rises above 7.5. In BEL Agr Bio however, liming opens this window for some time in each 10-year liming period because pH temporarily rises into the range 5–7. The phenomenon of renewed clay transport caused by agriculture has often been observed in loess soils (e.g. in the form of agricutans; Jongerius, 1970), although this is partly triggered by tillage practices and the occurrence of bare soil surface as well (Macphail et al., 1990). The leaching climate in the BEL scenarios caused a window for possible clay migration to open in the Bölling/Alleröd and to close (at least within 120 cm depth) at the start of the Holocene. The much dryer climate in the HUN scenarios causes the window to open at the start of the Atlanticum period. In the HUN scenarios a weak possibility for clay migration was calculated for the Bölling/Alleröd period. This may be connected to the flushing of soil when it is thawing. The scenario with absence of bioturbation shows a much larger time-window for possible clay migration. When no calcite is brought up, pH in the

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

475

Fig. 13. Simulated bulk density.

topsoil drops more rapidly thus opening the window for possible clay migration in the Alleröd. Soil morphological descriptions often yield important information on the genesis of a soil profile, but there exists little experience with modelling the formation of soil structure as a function of cryo- and bioturbation processes. Therefore an indicator for the preservation of frost-related soil morphology IPFMy was defined. This indicator is based on the assumption that if a given mass fraction of the soil is dissolved/eluviated, precipitated/illuviated or turbated in a year y, the same fraction of the soil volume will lose its original morphological characteristics due to fading, masking or erasure respectively.

On the other hand, freezing of the soil may cause new frostrelated morphology (e.g., a platy structure in the topsoil) and IPFMy should then be reset. IPFMy is calculated for each year and soil compartment by:  IPFMy ¼

1 IPFMy1  fturby  fcalcy  fgypsy  fclayy

8 T b  2 -C 8 T z  2 -C

ð8Þ xwhere fturby = (1 − bioturbation), fcalc y = (1 − abs(calcitey − calcite y − 1 )), fgyps y = (1 − abs(gypsum y − gypsum y − 1 )), fclay y = ICM y,pH × dispersible y) and T y is the soil temperature

Fig. 14. Simulated exchangeable K.

476

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

Fig. 15. Clay dispersion indicator based on simulated pH. A value of 1 = likely, 0 = unlikely dispersion.

(°C) at 31–12 in year y. Bioturbation, calcite, gypsum and dispersible are all mass fractions. For dispersible y the value of 0.0005 was arbitrarily chosen, thus it was assumed that (yearly) maximally 0.05% of the mass of a soil compartment may be dispersed and transported. A value of IPFM y of 1 indicates that a soil morphology related to freezing is likely to be found in year y. The development of IPFM is depicted in Fig. 16. The BEL scenarios indicate that most frost-related morphology (indicator value of 0.54) is left at approximately 45 cm depth at present (year 0) due to the combined effects of decarbonization, clay transport and bioturbation over time. Field evidence shows the presence of consolidated periglacial Btx-horizons starting at comparable depths in natural soils. These Btx were according to Van Vliet-Lanoë et al. (1992) and Brahy et al. (2000) formed under near-permafrost conditions in the Younger Dryas. Although no permafrost was simulated for

the Younger Dryas, IPFM appears to reproduce the depth of maximal evidence for frost-related morphology. Disturbance in the BEL Agr scenario is slightly higher than in the BEL Nat scenario due to the higher possibility of clay transport in BEL Agr. The HUN scenarios show less disturbance in the subsoil. Surprisingly, the scenario without bioturbation shows the highest disturbance of frost-related morphology. This is because decalcification, formation of a calcic horizon and the possibility for clay migration are more pronounced in this scenario. The HUN Agr Bio scenario appears to conserve morphology best, because decarbonization, formation of a calcic horizon and clay migration are less pronounced due to the liming. It should be noted that the indicator is not fieldverified and may be subject to error since the bioturbation history may be more variable in intensity and with depth than was input to the model and actual clay migration is not simulated. However, defining and using morphological

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

477

Fig. 16. Indicator for disturbance of frost-related morphology based on bioturbation, depletion and accumulation of calcite, the indicator for clay dispersion and the occurrence of frozen soil. A value of 1 = unlikely, 0 = likely disturbance.

indicators does open new possibilities for the verification of soil models.

are currently not included in the model and that may need to be implemented for a wider application domain:

4.3. Model performance and application domain issues

i The issue of not including the effect of strain in the model is assumed to be of limited importance in the studied scenarios since simulated differences in bulk density are not large (starting value 1.3, final range between 1.15 and 1.50). When clay migration is to be simulated, this situation may change; ii In the BEL scenarios, pH dropped to values near 4 and aluminium chemistry as well as weathering of soil minerals should be taken into account for its effects on pH and

Each scenario run took approximately three 100% CPU-days (six days runtime) on a Pentium 4, 3.4 GHz PC running under Windows XP, which is reasonable for small numbers of scenarios. All simulations were completed without software stability problems. Technical performance was therefore concluded to be good. There are a number of processes that

478

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479

occupation of the exchange complex (Hoosbeek and Bryant, 1994; De Vries et al., 1995); iii Simulation of clay transport will allow more realistic scenarios to be analysed, because the clay migration indicator does not predict the changes in soil porosity and associated soil physical characteristics that occur in Bt-horizons; iv Only part of the soil–vegetation interaction is simulated since vegetation is input to the model and is limited to 4 vegetation types; v The effect of erosion and sedimentation processes can currently be simulated only by adding or removing whole (5 cm) model compartments. To allow evaluation of more gradual effects, which is essential to reconstruct soilscapes, some flexibility in compartment thickness must be realized. 5. Conclusions The SoilGen1-model has been developed, applied and tested for the simulation of soil development in calcareous loess. Test results indicate that the model is of sufficient quality to be used as a deterministic temporal interpolator to reconstruct soil development in calcareous parent materials such as loess for a variety of climates. Verification of process models for soil development is highly laborious when only measured data at the current stage of development are available. Verification with process rates of long temporal extents such as decarbonization rates is therefore of great value. Incorporation of the bioturbation process in pedogenetic models is essential, as decarbonization in the topsoil is slowed down up to 1000s of years in continental climates and as there are indications that bioturbation affects the time-window for clay migration. The models' sensitivity to initial and time-dependent boundary conditions has only partly been explored by the scenarios, thus few statements can be done on issues like path way dependency or convergent soil development for various climate evolutions. To widen the application domain towards more acidic conditions, processes such as aluminium chemistry and mineral weathering should be added to the model. Useful other additions would be the simulation of clay transport and of soil–vegetation interactions. In case processes such as clay transport and development of soil structures cannot be simulated due to lack of process knowledge, indicators for the occurrence or magnitude can be obtained and compared to observations. These will give insight on how these processes can be incorporated into future models. Acknowledgements We are grateful to Dr. Basil Davis for supplying the January and July temperature data for central-western Europe and central-eastern Europe (Davis et al., 2003) and the precipitation data for the same regions (Davis, unpubl.), which were of great benefit. Thanks are also due to Dr. C. Verbruggen and Dr. E. van Ranst for providing references and data on the development of vegetation and loess soils in Belgium.

References Addiscott, T.M., Wagenet, R.J., 1985. Concepts of solute leaching in soils: a review of modelling approaches. Journal of Soil Science 36, 411–424. Brahy, V., Peeters, G., Iserentant, A., Deckers, J., Delvaux, B., 2000. Estimation of a soil weathering stage and acid neutralizing capacity in a toposequence Luvisol–Cambisol on loess under deciduous forest in Belgium. European Journal of Soil Science 51, 1–13. Brimhall, G.H., Dietrich, W.E., 1987. Consecutive mass balance relations between chemical composition, volume, porosity and strain in metasomatic hydrochemical systems: results on weathering and pedogenesis. Geochimica et Cosmochimica Acta 51, 567–587. Campbell, C.A., Cameron, D.R., Nicholaichuk, W., Davidson, H.R., 1977. Effects of fertilizer N and soil moisture on growth, N content and moisture use by spring wheat. Canadian Journal of Soil Science 57, 289–310. Coleman, K., Jenkinson, D.S., 2005. RothC-26.3: a model for the turnover of carbon in soil. Model Description and Users Guide. November 1999 Issue (Modified April 2005). http://www.rothamsted.bbsrc.ac.uk/aen/carbon/ download.htm (accessed January 2007). Dann, R.L., Close, M.E., Lee, R., Pang, L., 2006. Impact of data quality and model complexity on prediction of pesticide leaching. Journal of Environmental Quality 35 (2), 628–640. Davis, B.A.S., Brewer, S., Stenvenson, A.C., Guiot, J., Data Contributors, 2003. The temperature of Europe during the Holocene reconstructed from pollen data. Quaternary Science Reviews 22, 1701–1716. De Vries, W., Posch, M., 2003. Derivation of Cation Exchange Constants for Sand, Loess, Clay and Peat Soils on the Basis of Field Measurements in the Netherlands. Alterra-Report, vol. 701. Alterra Green World research, Wageningen. 50 pp. De Vries, W., Kros, J., Van Der Salm, C., 1995. Modeling the impact of acid deposition and nutrient cycling on forest soils. Ecological Modelling 79 (1–3), 231–254. DeNovio, N.M., Saiers, J.E., Ryan, J.N., 2004. Colloid movement in unsaturated porous media: recent advances and future directions. Vadose Zone Journal 3, 338–351. Dokuchaev, V.V., 1883. Russian Chernozem (Russkii Chernozem) (Translated from Russian by N. Kaner) Israel Prog. For Sci. Transl., Jerusalem, 1967. Egli, M., Fitze, P., 2001. Quantitative aspects of carbonate leaching of soils with differing ages and climates. Catena 46, 35–62. Gobat, J.-M., Aragno, M., et Matthey, W., 1998. Le sol vivant. Presses polytechniques et universitaires romandes, Lausanne. Hargreaves, G.H., Samani, Z.A., 1985. Reference crop evapotranspiration from temperature. Transaction of ASAE 1 (2), 96–99. Hoosbeek, M.R. and Bryant, R.B. 1994. Developing and adapting soil process submodels for use in the pedodynamic Orthod model. p. 111–128. In R.B. Bryant and R. Arnold (ed.) Quantitative modeling of soil forming processes. SSSA Spec. Publ. vol. 39. ASA, CSSA, and SSSA, Madison, WI. 185 pp. Hutson, J.L., 2003a. LEACHM Model description and user's guide. . http:// www.scieng.flinders.edu.au/cpes/people/hutson_j/leachweb.html (accessed January 2007). Hutson, J.L. 2003b. LEACHM — a process-based model of water and solute movement, transformations, plant uptake and chemical reactions in the unsaturated zone. Version 4. Research Series No R03-1, Dept of Crop and Soil Sciences, Cornell University, Ithaca, NY. Hutson, J.L., Wagenet, R.J., 1992. LEACHM: Leaching Estimation and Chemistry Model: A Process based Model of Water and Solute Movement Transformations, Plant Uptake and Chemical Reactions in the Unsaturated Zone. Continuum, vol. 2. Water Resources Inst., Cornell University, Ithaca, NY. Version 3. IUSS Working Group WRB, 2006. World Reference Base for Soil Resources 2006, 2nd edition. World Soil Resources Report, vol. 103. FAO, Rome. Jabro, J.D., Jabro, A.D., Fox, R.H., 2006. Accuracy and performance of three water quality models for simulating nitrate nitrogen losses under corn. Journal of Environmental Quality 35 (4), 1227–1236. Jalali, M., Rowell, D.L., 2003. The role of calcite and gypsum in the leaching of potassium in a sandy soil. Experimental Agriculture 39 (4), 379–394. Jenkinson, D.S., Coleman, K., 1994. Calculating the annual input of organic matter to soil from measurements of total organic carbon and radiocarbon. European Journal of Soil Science 45, 167–174.

P.A. Finke, J.L. Hutson / Geoderma 145 (2008) 462–479 Jenny, H., 1941. Factors of Soil Formation: A System of Quantitative Pedology. McGraw-Hill, New York. 281 pp. Jongerius, A., 1970. Some morphological aspects of regrouping phenomena of Dutch soils. Geoderma 4, 311–331. Kononova, M.M., 1975. Humus of virgin and cultivated soils. In: Gieseling, J.E. (Ed.), Soil components. I. Organic Components. Springer, Berlin, pp. 475–526. Kros, J., Pebesma, E.J., Reinds, G.J., Finke, P.A., 1999. Uncertainty assessment in modelling soil acidification at the European scale: a case study. Journal of Environmental Quality 28 (2), 366–377. Liu, X., Rolph, T., Bloemendal, J., Shaw, J., Liu, T., 1995. Quantitative estimates of palaeoprecipitation at Xifeng, in the Loess Plateau of China. Palaeogeography, Palaeoclimatology, Palaeoecology 113 (2–4), 243–248. Lundin, L.-C. 1989. Water and heat flows in frozen soils. Basic theory and operational modeling. Ph.D. thesis. Uppsala University, Uppsala, Sweden. Macphail, R.I., Courty, M.A., Gebhardt, A., 1990. Soil micromorphological evidence of early agriculture in north-west Europe. World Archaeology 22 (1), 53–69. McBratney, A.M., Mendonça Santos, M.L., Minasny, B., 2003. On digital soil mapping. Geoderma 117, 3–52. Moldrup, P., Olesen, T., Schjonning, P., Yamaguchi, T., Rolston, D.E., 2000. Predicting the gas diffusion coefficient in undisturbed soils fro soil water characteristics. Soil Science Society of America Journal 64, 94–100. Navrátil, T. 2003. Biogeochemistry of the II.A group elements in a forested catchment. PhD-thesis. Charles University, Prague. Parton, W.J., Scurlock, J.M.O., Ojima, D.S., Gilmanov, T.G., Scholes, R.J., Schimel, D.S., Kirchner, T., Menaut, J.-C., Seastedt, T., Garcia Moya, E., Kamnalrut, A., Kinyamario, J.I., 1993. Observations and modeling of biomass and soil organic matter dynamics for the grassland biome worldwide. Global Biogeochemical Cycles 7, 785–809. Petit, J.R., Jouzel, J., Raynaud, D., Barkov, N.I., Barnola, J.M., Basile, I., Bender, M., Chappellaz, J., Davis, J., Delaygue, G., Delmotte, M., Kotlyakov, V.M., Legrand, M., Lipenkov, V., Lorius, C., Pépin, L., Ritz, C., Saltzman, E., Stievenard, M., 1999. Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica. Nature 399, 429–436. Plinius the Young. Natural history (book 17, 44). Rawls, W.J., Brakensiek, D.L., 1985. Agricultural management effects on soil water retention. In: DeCoursey, D.G. (Ed.), Proceedings of the 1983 Natural Resources Modelling Symposium. U.S. Department of Agriculture, Agricultural Research Service, ARS-30. 532 pp. Rybnícknova, E., Rybnicek, K., 1996. Czech and Slovak Republik. In: Berglund, B.E., Birks, H.J.B., Ralska-Jasiewiczowa, M., Wright, H.E. (Eds.), Palaeoecological Events during the Last 15000 Years. Regional Syntheses of Palaeoecological Studies of Lakes and Mires in Europe. Wiley, Chistester.

479

Schaetzl, R.J., Schwenner, C., 2006. An application of the Runge “Energy Model” of soil development in Michigan's upper peninsula. Soil Science 171/2, 152–166. Seppä, H., Birks, H.J.B., 2001. July mean temperature and annual precipitation trends during the Holocene in the Fennoscandian tree-line area: pollen-based climate reconstructions. Holocene 11 (5), 527–539. Singer, A.C., Jury, W., Luepranchai, E., Yahng, C.-S., Crowley, D.E., 2001. Contribution of earthworms to PCB bioremediation. Soil Biology & Biochemistry 33, 765–776. Smith, P., Smith, J.U., Powlson, D.S., McGill, W.B., Arah, J.R.M., Chertov, O.G., Coleman, K., Franko, U., Frolking, S., Jenkinson, D.S., Jensen, L.S., Kelly, R.H., Klein-Gunnewiek, Komarov, A.S., Li, C., Molina, J.A.E., Mueller, T., Parton, W.J., Thornley, J.H.M., Whitmore, A.P., 1997. A comparison of the performance of nine soil organic matter models using datasets from seven long term experiments. Geoderma 81, 153–225. Stumm, W., Morgan, J.J., 1970. Aquatic Chemistry. Wiley-Interscience. Thompson, K., Parkinson, J.A., Band, S.R., 1997. A comparative study of leaf nutrient concentrations in regional herbaceous flora. New Phytologist 136, 679–689. Tillotson, W.R., Robbins, C.W., Wagenet, R.J., Hanks, R.J., 1980. Soil water, solute, and plant growth simulation. Bulletin 502 Utah State Agr. Exp. Stn., Logan, Utah. 53 pp. Van Breemen, N., Buurman, P., 2002. Soil Formation. Kluwer Acad. Publishers, Dordrecht. second, revised edition. Van Ranst, E., 1981. Genesis and properties of silty forest soils in central Belgium and the Ardennes (in Dutch). Unpublished PhD-thesis. Van Vliet-Lanoë, B., Fagnart, J.P., Langohr, R., Munaut, A., 1992. Importance de la succession des phases écologiques anciennes et actuelles dans la différentiation des sols lessivés de la couverture loessique d'Europe occidentale: argumentation stratigraphique et archéologique. Science du Soil 30, 75–93. Verbruggen, C., Denys, L., Kiden, P., 1996. Belgium. In: Berglund, B.E., Birks, H.J.B., Ralska-Jasiewiczowa, M., Wright, H.E. (Eds.), Palaeoecological Events during the Last 15000 Years. Regional Syntheses of Palaeoecological Studies of Lakes and Mires in Europe. Wiley, Chistester. Wierenga, P.J., Nielsen, D.R., Hagan, R.M., 1969. Thermal properties of a soil based upon field and laboratory measurements. Soil Science Society of America, Proceeding 33, 354–360. Wyszkowski, M., Wyszkowska, J., Włodkowska, L., 2006. Correlations between macroelements content of spring barley and the enzymatic activity of soil contaminated with copper, zinc, tin and barium. Electronic Journal of Polish Agricultural Universities 9 (2). http://www.ejpau.media.pl/volume9/ issue2/art-02.html.