SHIMMER (1.0): a novel mathematical model for microbial and biogeochemical dynamics in glacier forefield ecosystems

. SHIMMER (Soil biogeocHemIcal Model for Microbial Ecosystem Response) is a new numerical modelling framework designed to simulate microbial dynamics and bio-geochemical cycling during initial ecosystem development in glacier foreﬁeld soils. However, it is also transferable to other extreme ecosystem types (such as desert soils or the surface of glaciers). The rationale for model development arises from decades of empirical observations in glacier foreﬁelds, and enables a quantitative and process focussed approach. Here, we provide a detailed description of SHIMMER, test its performance in two case study foreﬁelds: the Damma Glacier (Switzerland) and the Athabasca Glacier (Canada) and analyse sensitivity to identify the most sensitive and un-constrained model parameters. Results show that the accumulation of microbial biomass is highly dependent on variation in microbial growth and death rate constants, Q 10 values


Introduction
Ice fronts in polar and alpine regions are retreating as a result of climate warming, and as a consequence, glacier forefield areas in high-latitude and high-altitude regions are rapidly expanding (Graversen et al., 2008;ACIA, 2005).Glacier coverage in upland Alpine regions in Europe has declined by up to 30 % from the 1970s to 2003, exposing roughly 860 km 2 of previously ice-covered land area (Paul et al., 2011).Similarly, rapid glacial retreat has been observed in Iceland (Staines et al., 2014), North America (Insam and Haselwandter, 1989;Hahn and Quideau, 2013;Ohtonen et al., 1999;Sattin et al., 2009), Asia (Liu et al., 2012) and Svalbard (Moreau et al., 2008).These vast expanses of newly uncovered land have been locked under ice for tens of thousands of years, are typically highly oligotrophic, with low nutrient budgets and are subject to harsh and rapidly fluctuating environmental conditions.They potentially play a significant yet largely unexplored role in large-scale biogeochemical cycling and climate (Dessert et al., 2003;Anderson et al., 2000;Smittenberg et al., 2012;Berner et al., 1983), global methane budgets (Kirschke et al., 2013), the global phosphorus cycle (Filippelli, 2002;Follmi et al., 2009) and the productivity of downstream and coastal ecosystems (Anesio et al., 2009;Mindl et al., 2007;Fountain et al., 2008;Anderson et al., 2000).Furthermore, the initial stages of Published by Copernicus Publications on behalf of the European Geosciences Union.microbial community development in soils are fundamental to understanding life in extreme environments, which may serve as an analogue for survival and habitability of environments currently assumed devoid of life on this planet and others.
Microbial communities are the primary colonisers of recently exposed soils, and are thought to be fundamental in soil stabilisation, the build-up of carbon and nutrient pools, and facilitating the establishment of higher-order plants (Schulz et al., 2013;Bradley et al., 2014).A conceptual overview of forefield nutrient cycling is presented in Fig. 1.Recently de-glaciated soils vary in their mineralogical and microbial compositions.Total organic carbon (TOC), total nitrogen (TN), and total phosphorus (TP) content of newly exposed glacial forefield soils is low, in the range of 0.1-40.0mg C g −1 , 0.1-2.0mg N g −1 , and 2-8 µg P g −1 (Bradley et al., 2014).However, these concentrations typically increase by 1 to 2 orders of magnitude with ageing from newly exposed to well-developed (decades old) soil.Three distinct sources contribute carbon to recently exposed soils: autochthonous primary production by autotrophic microorganisms, allochthonous material deposited on the soil surface (from wind, hydrology, biology and precipitation) and ancient organic pools derived from under the glacier.Organic matter accumulation from all three sources supports the development of heterotrophic communities, yet their relative significance remains unknown.The continual autotrophic production, heterotrophic re-working and allochthonous deposition lead to the accumulation of organic material, which supports higher abundances and diversity of microorganisms (Bradley et al., 2014;Schulz et al., 2013).Nitrogen is derived from active nitrogen-fixing organisms, allochthonous deposition and degradation of organic substrate.Bioavailable phosphorus is usually abundant in the topsoil or bedrock of glaciated regions from weathering of the mineral surface, and can also be liberated due to the degradation of organic molecules.
Hypotheses relating accumulation of carbon and nutrients to increasing species richness and diversification have tended to be descriptive and qualitative, rather than quantitative.Little is known about the main drivers of microbial succession and controls on abundance, diversity and activity, and this limits our understanding of glacier forefield contribution to global nutrient pathways and our ability to predict how these rapidly expanding ecosystems may respond in the future, including their potential impact on atmospheric CO 2 , global climate and downstream productivity.This lack of understanding can be partly attributed to the difficulty of quantifying the different external organic matter and nutrient fluxes, as well as disentangling the complexity of biogeochemical processes underlying the microbial dynamics and soil carbon and nutrient build-up along the chronosequence by observations and/or laboratory experiments alone.Numerical models are useful tools in this context as they can not only help to disentangle the complex process interplay, diagnose the fluxes between ecosystem components and, ultimately, predict the sensitivity and response of an environment to changing environmental and climatic conditions, but also help identify important data and knowledge gaps and hence guide the design of efficient field campaigns and laboratory studies directly targeted at closing these gaps.Nevertheless, a modelling framework that could be used to explore microbial dynamics and associated nutrient cycling in glacier forefields currently does not exist.
The development of soil models has been common in the past and important in informing soil management, policy and prediction (McGill, 2007(McGill, , 1996)), for example in understanding the contribution of soil organic matter (SOM) to the formation of stable aggregate soils, the ease of soil cultivation, water holding characteristics and the risk of physical damage and compaction.The explicit inclusion of soil microbial dynamics has been shown to drastically improve the performance of these models (Wieder et al., 2013).There are many different types of soil models in use today across a range of scales and purposes, such as informing agricultural policy, understanding biogeochemical cycling and soil food webs and the feedbacks between soil processes, hydrology and the atmosphere (Stapleton et al., 2005;Blagodatsky and Richter, 1998;Knapp et al., 1983;Grant et al., 1993;German et al., 2012;Ingwersen et al., 2008;Leffelaar and Wessel, 1988;Kuijper et al., 2005;Kravchenko et al., 2004;Parton et al., 1988;Garnier et al., 2001;Darrah, 1991;Foereid and Yearsley, 2004;Vandewerf and Verstraete, 1987b;Long and Or, 2005;Maggi and Porporato, 2007;Moorhead and Sinsabaugh, 2006;Panikov and Sizova, 1996;Toal et al., 2000;Zelenev et al., 2000;Scott et al., 1995).However, although these models include an explicit microbial component, SOM models are tailored towards research questions that are focussed on geochemistry and specifically organic matter dynamics rather than biology.Forefield ecosystems are characterised by extreme and highly variable environmental conditions and rapidly changing compositions of microbial communities whose interplay results in unique chronosequence dynamics (Bradley et al., 2014).There is not a single model that can represent the unique forefield development without an unacceptable level of abstraction and simplification of the system.
Therefore, we developed the new model framework SHIMMER (Soil biogeocHemIcal Model for Microbial Ecosystem Response) to quantitatively simulate the initial stages of ecosystem development and assess biogeochemical processes in the forefield of glaciers.The code is written and executed in the free open-source computing environment and programming language R, which is available to download on the web (http://www.r-project.org/).The current version of the model is designed to represent the microbial community prior to the establishment of plants, and therefore only the initial stages of chronosequence development will be assessed.Microbial communities may be heavily structured by establishing vegetation (Brown and Jump-Figure 1.A conceptual overview of the main transformations and fluxes of nutrients and organic matter in and along a typical de-glaciated forefield.
Table 1.Outline of quantitative and predictive modelling strategy.

Quantitative modelling Details
Do initial microbial communities rely on autochthonous or allochthonous carbon?
Explore scenarios of carbon loading to determine the reliance on allochthonous carbon vs. autochthonous production.Assess the role of nitrogen as a key limiting nutrient.
Quantify nitrogen budget to assess relative importance of nitrogen fixation vs. DIN assimilation in recently exposed and old soils.Assess how microbial diversity affects soil development.
Test various community assemblages and characterise the soil environment they create.

Predictive modelling Details
Identifying tipping points within the system.Do population dynamics, external nutrient loading, seasonal changes, disturbance events and climate change trigger significant changes in soil development?Assess importance of seasonality.
Explore the dynamics of the "non-growing season", where there is little observational data due to inaccessibility and difficult field work conditions.Assess importance of disturbance events and community reset.
Test the effect of hydrological disturbance in proglacial zone (rich in nutrients and organic matter, but may reset established communities and enhance leaching of substrate) Assess sensitivity to future climate change Explore the effect of climate and anthropogenic impact using scenario-based predictions.Draw attention to gaps in our understanding and areas of future research.Models are driven primarily by empirical relationships and observational data; henceforth it is likely to become apparent where future fieldwork and lab work efforts should be focussed.ponen, 2014), and the physical properties of vegetated soils are considerably different in terms of water retention, ultraviolet exposure, temperature fluctuations (Ensign et al., 2006;King et al., 2008) and nutrient status (Kastovska et al., 2005;Schutte et al., 2009).
This paper provides a comprehensive description of the modelling framework.A first model performance test is conducted on the basis of existing published data from the Damma Glacier forefield in Switzerland (Bernasconi et al., 2011;Brankatschk et al., 2011;Guelland et al., 2013b) and from the Athabasca Glacier in Canada (Insam and Haselwandter, 1989).The newly developed model is then used to conduct an extensive parameter sensitivity study.SHIM-MER is not only a new model framework but also part of an interdisciplinary, iterative, open-source, model-data-based approach fully integrating fieldwork and laboratory experiments with model development, testing and application.The model scope and complexity of the first version of SHIM-MER is informed by a comprehensive review of glacier forefield research (Bradley et al., 2014) and the data set collected during a first field campaign to characterise the forefield of a retreating Svalbard glacier.The model is kept as general as possible, and thus is easily transferable to other microbial ecosystems such as desert soils, ice surfaces (e.g.cryoconite), microbial mats and the built environment (e.g.fuel and chemical storage).The model is dynamically sufficient; i.e. the minimum processes that are needed to resolve the system and provide useful output are included.In addition, it is intended that new model developments will guide and inform future field and laboratory studies so that subsequent versions of the model will run with narrower plausible ranges of parameters, and explicitly resolve processes that currently cannot be constrained on the basis of available data.SHIM-MER can thus be considered as the first step of an interdisciplinary, iterative approach (illustrated in Fig. 2 and Table 1) where new data inform new model developments that will result in new insights, which in turn will inform new laboratory and field experiments, etc.It thus not only contributes to more accurate quantitative predictions that enable a deeper understanding of the processes that control microbial communities, their role on global biogeochemical cycles and their response to climate variations, but also provides an ideal platform for the synthesis and exchange of knowledge and information across different disciplines.
The final developed model presented here is: -structurally (i.e.spatial resolution, number of species, processes included) and mechanistically (i.e.process formulation) complex enough to describe the required properties for carbon, nitrogen and phosphorus turnover necessary to address the questions identified in the main aims (Table 1); -structurally and mechanistically simple enough to constrain and validate parameters and simulation results by available data and literature; -able to operate with numerical efficiency on various timescales from days to decades, in order to represent an entire chronosequence development in sufficient detail; -applicable to a range of environments; -structurally stable, conserves mass and provides robust numerical output.

Model description
The following sections provide a detailed description of SHIMMER and its implementation.A conceptual model diagram is presented in Fig. 3.The state variables are listed in Table 2. Derived variables are included in the model to quantify mass budgets and transfers between pools.These are listed in Table 3, along with their formulation.The mathematical formulation of the model is presented in Table 4. Parameters are summarised in Table 5.

Physical support
The model support of SHIMMER (1.0) is 0-D; i.e. there is no specific spatial discretisation (e.g.depth).This choice is informed by the quality of observational data.

Model components
The model contains pools of microbial biomass, organic matter and both dissolved inorganic and organic nitrogen and phosphorus (Table 2).A system of coupled ordinary differential equations describes the transfers and transformations of these pools (Table 4).

Microbial dynamics
In its current form, SHIMMER distinguishes between autotrophic (A) and heterotrophic (H ) microbial communities, which are further subdivided into three categories (A 1−3 , DIN consumed by all microbes DIN Consumed cum_n f DIN fixed by A 3 and referred to as "soil microbes".A 3 and H 3 are "nitrogenfixing microbes" or "nitrogen fixers", able to fix atmospheric N 2 gas as a source of nitrogen when dissolved inorganic nitrogen (DIN) stocks become limiting, giving them a competitive advantage.Examples of A 3 and H 3 are Nostoc and Rhizobiales, respectively.The overall rate of change of autotrophic and heterotrophic communities is described by Microbial dynamics are thus governed by three terms: growth (U ), losses (G), and the rate of exudate and exopolymeric substances (EPS) production (X).

Growth U
Microbial biomass is created by chemolithoautotrophic (U A1 ), photo-autotrophic (U A2 and U A3 ) and heterotrophic (U H 1−3 ) growth.Subglacial chemolithoautotrophs (A 1 ) acquire their energy through chemical synthesis of mineral compounds.Autotrophic soil species (A 2 and A 3 ) acquire energy through photosynthesis, and accordingly require light in order to grow.Chemolithoautotrophic growth (Eq.3) and photo-autotrophic growth (Eq.4) are described as where autotrophic growth is calculated according to a maximum growth rate (I max ) (at temperature (T ) = T ref ) which is further affected by substrate availability, temperature, photosynthetically available radiation (PAR), and nutrient concentrations (Mur et al., 1999;Van Liere and Walsby, 1982).Temperature dependency is described by a temperature response factor (T f ) with a Q 10 formulation (Table 4) effectively slowing down or speeding up all life processes (Soetaert and Herman, 2009;Yoshitake et al., 2010;Schipper et al., 2014) Rate of change of heterotrophic biomass (H 1−3 ) Rate of change of labile carbon substrate (S 1 ) Rate of change of refractory carbon substrate (S 2 ) Rate of change of labile organic nitrogen (ON 1 ) Rate of change of refractory organic nitrogen (ON 2 ) Rate of change of labile organic phosphorus (OP 1 ) Rate of change of refractory organic phosphorus (OP 2 ) Growth of N-fixing soil autotrophs (A 3 ) with nitrogen fixation Growth of N-fixing soil autotrophs (A 3 ) with DIN Growth of subglacial heterotrophs (H 1 ) from labile substrate Growth of soil heterotrophs (H 2 ) from labile substrate Growth of N-fixing soil heterotrophs (H 3 ) from labile substrate and nitrogen fixation Growth of N-fixing soil heterotrophs (H 3 ) from labile substrate and DIN Growth of subglacial heterotrophs (H 1 ) from refractory substrate Growth of N-fixing soil heterotrophs (H 3 ) from refractory substrate and nitrogen fixation Growth of N-fixing soil heterotrophs (H 3 ) from refractory substrate and DIN Growth of N-fixing soil heterotrophs (H 3 ) from labile substrate Growth of N-fixing soil heterotrophs (H 3 ) from refractory substrate

Rate of exudate and EPS production by autotrophs
Rate of exudate and EPS production by heterotrophs

Phosphorus cycle
Labile organic phosphorus degraded • (S 2Consumed ) Labile organic phosphorus accumulated OP 1Accumulation = P C • (G L + X T ) Refractory organic phosphorus accumulated Environmental and scaling equations Temperature factor response species (A 1 and H 1 ) are represented in SHIMMER as more energy conserving and adapted to harsh environmental conditions.Accordingly, their maximum growth rate (I max ) is reduced (by a factor p Sub ), and their nutrient limitation response (expressed by the half-saturation constants K S , K N and K P ) is reduced (by a factor k Sub ).

Y H
Growth efficiency of heterotrophs g C (g C consumed) −1 0.200 (Scott et al., 1995) 0.134 (German et al., 2012) 0.848 (Blagodatsky and Richter, 1998) --  kinetics (Table 4) (Holl and Montoya, 2005;Rabouille et al., 2006;Goebel et al., 2008).To account for the additional energy expenditure of nitrogen fixation, the growth efficiency (Y A and Y H ) and maximum growth rates (I max ) of nitrogen fixers whilst fixing nitrogen (U A3_N2 and U H 3_N2 ) are reduced by a factor n f (Breitbarth et al., 2008;LaRoche and Breitbarth, 2005).Whilst nitrogen fixers are actively fixing atmospheric nitrogen, soil DIN concentration is not a limiting factor on their growth rate and accordingly, the growthlimiting half-saturation expression for soil DIN (K n ) is discounted from U A3_N2 and U H 3_N2 .Similarly, whilst nitrogen fixers are using N 2 gas as their source of nitrogen, there is no uptake of DIN from the soil DIN pool (Table 4).Nitrogen fixation in the SHIMMER model is sensitive to many of the environmental factors which are often cited, including surrounding DIN concentrations, temperature and carbon and phosphorus limitation (Liu et al., 2011) (see e.g.U A3_N2 in Table 4).Heterotrophic microbes (H 1−3 ) acquire their energy through the degradation of organic substrate (S 1 and S 2 ) rather than by photosynthetic or chemolithoautotrophic processes.Their growth from labile (U H 2L ) and refractory (U H 2R ) organic matter is formulated in a similar way to Eqs. (3) and (4) but depends on the bioavailability of organic matter rather than light: The model distinguishes between two pools of organic matter.A reactive pool (S 1 ) comprises highly available and fresh organic compounds that are preferentially degraded by microorganisms and therefore may respond quickly to changing external conditions and inputs.A less reactive pool (S 2 ) represents the bulk of substrate present in the non-living organic component of soil.This organic matter degrades over longer timescales and therefore accumulates and may respond more slowly to changes in the environment.In order to express a preference of labile substrate, the parameters J S1 and J S2 (with J S1 > J S2 ) represent factors that scale the maximum rate at which labile carbon substrate (S 1 ) and refractory substrate (S 2 ) are utilised, respectively.In natural environments, most microorganisms live under fluctuating conditions, and as such their growth is inhibited in response to suboptimal conditions (Cowan et al., 2004).Organisms commonly reduce their metabolic activity and lower their energy expenditure in order to endure adverse environmental conditions.Accordingly, a large fraction of microorganisms in any environmental sample are in a metabolically inactive state (Lennon and Jones, 2011).The active fraction of microbial biomass is represented by the parameter d.This parameter is affixed to the growth and loss expressions.

Loss G
The loss terms (G Ai and G H i ) represent the net removal of biomass from the living biomass pools (A 1−3 and H 1−3 ).These are integral measures of natural death, viral lysis and grazing, which are lumped together in order to reduce complexity, to keep the number of parameters in a feasible and manageable range and appropriate to the current availability of data and understanding of the system.Mortality due to predation is usually density dependent (Kaitala et al., 1999;Levin, 1998).Accordingly, loss terms (G Ai and G H i ) are density dependent and are also sensitive to variations in soil temperature (T ), the active fraction d and adjustable rate parameters α A (autotrophs) and α H (heterotrophs): It is assumed that losses (G Ai and G H i ) form insoluble microbial necromass (organic matter), comprising of organic Geosci.Model Dev., 8, 3441-3470, 2015 www.geosci-model-dev.net/8/3441/2015/carbon substrate (S), organic nitrogen (ON) and organic phosphorus (OP), which enters the surrounding soil and is immediately available as substrate for heterotrophic growth.

Exudate production X
Exudate production by autotrophs (X Ai ) and heterotrophs (X H i ) (Eqs. 9 and 10) is proportional to growth rates (Allison, 2005), and is modulated by the exudation constants ex A (autotrophs) and ex H (heterotrophs): Exudate is highly reactive; therefore these stocks directly contribute to the labile substrate pools (S 1 , ON 1 , OP 1 ).

Organic matter dynamics
Organic matter dynamics are described by the following equations: Organic matter accumulates in the soil due to microbial loss (G Ai and G H i ), exudate and EPS production (X Ai and X H i ), and allochthonous external carbon substrate inputs (I S1 and I S2 ), and is depleted due to heterotrophic growth (U S1H i and U S2H i ) and leaching (W Si ).Substrate leaching (W Si ) is proportional to the mass of substrate and an adjustable parameter g Si .The coupling of substrate degradation to biomass growth (U S1H i and U S2H i ) is governed by the yield Y H , describing bacterial growth efficiency (biomass increase per unit of carbon substrate consumed).The remainder is respired as dissolved inorganic carbon (DIC).The adjustable parameter q represents the partitioning of substrate into the labile fraction (q) and refractory fraction (1 − q).

Nitrogen and phosphorus
The model accounts for DIN and DIP, as well as ON and OP species.The nitrogen and phosphorus components of living and dead microbial biomass are stoichiometrically coupled to microbial carbon (A 1−3 and H 1−3 ) and organic matter (S 1 and S 2 ) pools, respectively, by the N : C and P : C ratios of soil bacteria in Table 5 (Bernasconi et al., 2011).DIN dynamics are described as DIN is consumed by microbial growth and recycled through the heterotrophic degradation of organic matter.Atmospheric DIN deposition (I DIN ) and leaching (W DIN ) are potentially important processes in the nitrogen cycle in forefield soils (Williams et al., 2009;Brooks and Williams, 1999;Schimel et al., 2004).Phosphorus dynamics are represented in an almost identical way (however, note the absence of atmospheric fixation).

Numerical solution
The model is written and solved in the R statistical environment, an open-source computing environment and programming language.Due to the non-linear and complex nature of the equations which comprise the model, they must be solved numerically rather than analytically.SHIMMER uses the adaptive time-step solver "lsoda" from the deSolve package (Soetaert et al., 2010) to calculate the numerical solution.
Results are provided at daily time steps.On a standard desktop computer running R, the model usually takes less than 1 min to simulate 10 years of succession.

Forcing and initial conditions
The following external forcings drive and regulate the system's dynamics: -PAR (wavelength of approximately 400-700 nm) -nutrient leaching rate.
The presence of snow attenuates sunlight and inhibits PAR from reaching the soil surface.This is accounted for in preprocessing of forcing data according to the equation: whereby n is the irradiance (W m −2 ), x is the snow depth (m) and m is the extinction coefficient for snow (m −1 ).The extinction coefficients for various types of snow can be measured and an estimate of 6 m −1 is used in this instance to represent snow in glacier forefields (Greenfell and Maykut, 1977).External forcing and initial conditions for the Damma Glacier and Athabasca Glacier are presented in Fig. 4 and   (black).Forcing data are provided by the WSL (Switzerland), Achuff and Coen (1980) and Staveley AAFC (Canada).
4 Model runs

Nominal
The model is run with nominal parameters for a period representing 75 years of succession (the approximate length of the test data sets), starting on the 1 January, in order to provide a baseline model output from which parameters are varied to determine sensitivity.Leap years are ignored.The model is forced with meteorological data collected from the Damma Glacier, Switzerland (case study 1 -see Sect.4.3.1 for details).

Sensitivity and uncertainty analysis
A sensitivity study involving 24 model parameters is carried out to assess the stability of model output given variation in all key parameters.This is important when accounting for uncertainty, since high sensitivity of key parameters that have a relatively wide plausible range of values would lead to large uncertainties in predictions.Sensitivity analysis is considered across all state variables to assess the extent to which parameter variation influences whole model behaviour or only single variables.The model is run for 75 years under each sensitivity scenario, starting on 1 January, forced with data from the Damma Glacier (since results can be interpreted alongside the more detailed contextual observations).
The following model output (X) is explored: -total autotrophic biomass (average over the entire simulation run) -total heterotrophic biomass (average over the entire simulation run) Plausible ranges for model parameter values are constrained from values in published literature (Table 5).Many of the parameters show considerable variation, but the most confident values (their applicability to the glacier forefield system, the method by which the value was determined, and their occurrence in the literature) are used as nominal values.To explore sensitivity, uncertainty and linearity, plausible ranges are split into tenths, and simulations are run sequentially through all eleven possible values for each parameter.
Sensitivity around nominal values is quantified using a variation on the method presented in Xenakis et al. (2008).The relative sensitivity (λ) of a certain model output (X) to a parameter (p) is estimated according to where p is the nominal value for the parameter, X is the model output from nominal parameter values, δp is the difference in parameter value either side of the nominal value, and δX is the change in model output (simulated over the range of parameters identified in δp).For clarity, this is illustrated in Fig. 5a.The sensitivity (λ) quantifies the relation between the model output and variation in a single parameter as a first derivative of their relationship either side of the nominal value, and is normalised based on the magnitude of parameter and model output values.Thus, λ indicates the sensitivity of model output to parameter variation and also the direction (sign) of the change.A positive λ value indicates that an increase in the parameter value yields an increased value in the model output, whereas a negative λ value indicates that an increase in the parameter causes a decrease in the value of the model output.Values of λ further from zero indicate that the model output is highly sensitive to variation in the parameter.Model output is assessed graphically for each parameter (e.g.Fig. 6).First, the shape of the model output variation is assessed to see if the value for λ is representative of sensitivity.An unrealistic λ may be calculated if the nominal parameter is near a vertex and the variation in model output either side of the nominal value has an opposite sign (i.e. a parabola).This is illustrated in Fig. 5b, whereby δX is low, and thus a low λ value is obtained, even though the sensitivity is relatively high (i.e.X depends strongly on p).Second, each plot is assessed and a linear (e.g.Fig. 6a) or non-linear (e.g.Fig. 6b) relationship is attributed to each parameter.Finally, non-linear results are assessed to determine if the high- est sensitivity is around the nominal parameter value, since this has implications in interpreting the model output.If parameters are most sensitive near to the nominal values, there is a higher potential variation in model output and therefore potentially greater uncertainty in interpreting results.
To explore uncertainty (ø) associated with each parameter, the percentage variation in model output is calculated according to www.geosci-model-dev.net/8/3441/2015/Geosci.Model Dev., 8, 3441-3470, 2015 where X max and X min are the highest and lowest values for model output over the entire plausible range in parameter variation, and X is the model output with nominal parameters.

Optimisation
Model parameters implicitly account for all processes that are not explicitly accounted for in the model and, therefore, may vary across different environments.Based on the outcome of the sensitivity and uncertainty analysis, parameters are adjusted in an optimisation exercise to improve model fit to the validation data sets.Parameters are varied incrementally to determine the effect on the accumulation of microbial biomass and mean squared error is calculated with each model run.Parameters J S1 and J S2 (the relative bioavailability of labile and refractory substrate) are artefacts of the SHIMMER modelling structure.Therefore, these two parameters are the primary free parameters, which are adjusted to reduce mean squared error.Once a known optimum range for these parameters has been determined, the parameters that bear the highest sensitivity, uncertainty and non-linearity are adjusted.Parameters for which there is a high degree of confidence (narrow plausible ranges, lower sensitivity, linear behaviour and low uncertainty) are not adjusted in the optimisation exercise, since even relatively large changes in their value would cause only a small change in model output.
Given the wealth of physical, biological, genomic and chemical data available for the Damma Glacier, the focus of the analysis of model dynamics is on this data set.However, data from the Athabasca Glacier forefield provide additional support that the model can respond dynamically to predict the development of soils from a range of environments and study sites.This data set is more representative of the quality of data that is typically available for de-glaciated forefields.
The forefield chronosequence is roughly 650 m long and represents a range of soil ages from 0 years old to around 120 years old (Brankatschk et al., 2011).The underlying bedrock is mainly granitic gneiss (Frey et al., 2010) with a silty, sandy soil texture (Lazzaro et al., 2009).The site has a northeast exposition (Bernasconi et al., 2011) and an inclination of 25 % (Sigler and Zeyer, 2002).Soil pH decreases from pH 5.1 in initial soils (10 years ice free) to 4.1 in developed soils (ice free for 2000 years) and water holding capacity increases from 26 to 33 % (Brankatschk et al., 2011).Recently exposed sites at the Damma Glacier (ice free for 6 to 13 years) are characterised by mostly unvegetated, sandy-silty sediment, gravel and large rocks.Intermediate soils (ice free for 60 to 80 years) are characterised by increasing vegetation cover and soil structure resembling a typical soil profile.The old sites (ice free for roughly 120 years) are fully vegetated, with clearly visible soil horizons (Bernasconi et al., 2011).Molecular characterisation suggests that both specialised heterotrophs (α-, β-, γ -Proteobacteria), autotrophs (Cyanobacteria) and other nitrogen-fixing microbes are found in all samples from all ages (Duc et al., 2009).There is a clear increase in TOC with soil age (Bernasconi et al., 2011) from around 700 µg C g −1 in recently exposed soils to around 30 000 µg C g −1 in developed soils.Similarly, microbial biomass, TN and phosphorus increase by roughly an order of magnitude from recently exposed soils to developed soils (Bernasconi et al., 2011;Goransson et al., 2011).
The model is evaluated using least-squares error against four chemical analyses presented in the BigLink data set (Bernasconi et al., 2011): presented as C mic ; -carbon substrate (S 1 + S 2 ): calculated as TOC-C mic ; -ON (ON 1 + ON 2 ): calculated as TN-N mic ; -available DIP: presented as P resin .
Observational data were collected on 7 September (day 244 of the year), and are therefore compared to model output from day 244 of each year.The omission of sites older than 77 years (due to vegetation influence) leaves 16 samples ranging from 5 years to 77 years since ice retreat.The 5year data are used as initial conditions, leaving 3 data points in the "early soils" category and 12 from later stages of succession where there is relatively high plant abundance.Leastsquare error calculation and minimisation of errors are done only on those data points.(Goransson et al., 2011).Initial microbial biomass is assumed to be evenly distributed between all microbial groups of autotrophs and heterotrophs, and initial substrate bioavailability is assumed to be 40 % labile and 60 % refractory.Initial values for OP were not presented in the BigLink data set (Bernasconi et al., 2011), but were assumed to follow a stoichiometric ratio (Bernasconi et al., 2011).An initial value for DIN was taken from Brankatschk et al. (2011).
PAR, snow depth and soil temperature at 3 cm depth (collected by an automatic weather station in the Damma Glacier forefield) were provided by the WSL Institute for Snow and Avalanche Research SLF, Switzerland.Light intensity is provided in units of photons (µmol m −2 s −1 ) which are converted to PAR (W m −2 ) by a conversion factor (0.219).PAR (W m −2 ), snow depth (m) and soil temperature ( • C) are averaged to provide daily forcing to the model, and linear interpolation is used between any (very infrequent) missing data points.The seasonal data set is repeated for the duration of the model run (75 years) (Fig. 4).
Inputs of OP are assumed to be stoichiometrically linked to carbon substrate according to the measured C : P ratio of microbial biomass (Bernasconi et al., 2011).Allochthonous substrate input is assumed to be 30 % labile and 70 % refractory.Several additional assumptions are required to convert units of deposition per surface area to units of weight.When considering a 1cm deep soil profile, 1 g dry soil occupies a surface area of 0.869 cm 2 (Guelland et al., 2013b).Since substrate material and DIN is liberated when the snowpack melts, the yearly accumulation is prescribed evenly over 10 days when there is significant snowmelt: days 158-167 (7-16 May).DIP is typically liberated by rock weathering; however Frey et al. (2010) analysed the minerals liberated from the weathering of the granitic Damma Glacier bedrock material and did not find any traceable amounts of phosphorus.Different mineralogy is likely to considerably alter the importance of rock weathering as a source of phosphorus between locations, increasing the uncertainty for the amount of DIP generated by weathering processes.The annual input of DIP is prescribed as 80 µg P cm −2 yr −1 (equal to DIN input), but this release is spread evenly over the first snow-free months of each year, from day 167 to 206 (16 June-25 July).Prescribed allochthonous inputs are presented in Table 6.The proportion of the allochthonous nutrient input that is available to the soil represented by the model is adjusted by pa-

Case study 2: Athabasca Glacier, Canada
Published data from the Athabasca Glacier forefield, Canada (52.2 • N, 117.2 • W), are used as a second case study in the validation exercise (Insam and Haselwandter, 1989).The Athabasca Glacier forefield is a high-altitude (2740 m) site with soil ages from 5 to 225 years.The mineralogy is medium textured, mostly calcareous and neutral to slightly alkaline pH (Insam and Haselwandter, 1989).The Athabasca glacier forefield is less intensively studied, and accordingly there is less contextual information on the biogeochemical development of soils than the Damma Glacier.However, the soils in the earlier stages of development (< 100 years) provide a robust test of model behaviour and underlying system dynamics due to the sparseness of vegetation and lack of interference in the microbial signal from vascular plants.
The model is evaluated using least-squares error against two observed bulk biogeochemical variables (Insam and Haselwandter, 1989): Observational data were collected in July and are compared to model output from day 196 of each year.Sites older than 50 years should be interpreted cautiously due to the influence of establishing vegetation.The 5-year data are used as initial conditions.Initial microbial biomass is assumed to be evenly distributed between all microbial groups of autotrophs and heterotrophs, and initial substrate bioavailability is assumed to be 40 % labile and 60 % refractory.Since there are no quantitative estimates of DIN, DIP, ON and OP, initial inorganic nutrient concentrations are assumed to follow the same ratio as the Damma Glacier case study, and organic material follows a stoichiometric ratio (Bernasconi et al., 2011).Annual profiles of monthly average soil temperature (at 5 cm depth) and snow depth are obtained from published literature (Achuff and Coen, 1980) 4).
provide daily forcing data (Fig. 4).Daily solar irradiance data from 2014 are obtained from the Alberta Agriculture and Rural Development Agroclimatic Information Service for a nearby meteorological station (Stavely AAFC, 50.2 • N, 113.9 • NW; 1360 m).These are repeated year-on-year for the duration of the model run.There is no observational, experimental or modelled data of sufficient quality to provide forcings of allochthonous inputs to the Athabasca Glacier forefield.Therefore estimations from the Damma Glacier are used and parameters v Sub , v DIN and v DIP are adjusted.As with case study 1, optimisation is carried out based on the results of the sensitivity study, and a minimisation of leastsquared error and visual fit to data is carried out based on numerous model runs varying parameters that were identified in the sensitivity test.

Nominal
The model behaviour under nominal parameters, and forced with meteorological data and initial conditions from the Damma Glacier is presented in Fig. 7. Total microbial biomass is initially stable at roughly 3.4 µg C g −1 (year 1), followed by an exponential growth phase to 46.8 µg C g −1 (year 15), and then a decline to near-steady-state around 31.0 µg C g −1 , varying seasonally by roughly 12.0 µg C g −1 .
Geosci  By the final year of simulation (year 75), the community has evolved (from an even split 16.7 % per pool) such that the most dominant pool is soil autotrophs (A 2 = 35.8%), followed by nitrogen-fixing autotrophs (A 3 = 27.4 %) and subglacial chemolithoautotrophs (A 1 = 12.3 %), with all heterotrophic biomass (H 1−3 ) making up the remaining 24.5 %.Total bacterial production rises steadily from 0.3 µg C g −1 yr −1 (year 1) to 114.2 µg C g −1 yr −1 (year 15), after which (year 31 onwards) bacterial production declines by roughly a half.Autotrophs are consistently the highest producers, responsible for between 72.6 and 89.2 % of the total bacterial production.
There is a steady accumulation of carbon substrate throughout the entire simulation, from 735.4 µg C g −1 (year 1) to 4129.2 µg C g −1 (year 75); however substrate becomes more refractory (39.4 % labile in year 1, 3.7 % labile in year 75).ON and OP follow similar dynamics to carbon substrate (S 1 and S 2 ).The accumulation of substrate is derived from autotrophic activity, the build-up of microbial necromass and allochthonous deposition.DIN and DIP accumulate during the first 14 years of the simulation, after which DIN increases and DIP declines.

Sensitivity
Sensitivity analysis is presented in Fig. 8.The accumulation of autotrophic and heterotrophic biomass is most sensitive to variation in Q 10 (λ ≥ 0.70).Biomass accumulation is also highly sensitive to adjustments in the active fraction (d) (−0.55 ≤ λ ≤ −0.52), the bioavailability of refractory substrate (J S2 ) (0.44 ≤ λ ≤ 0.64), the partitioning of microbial necromass into labile and refractory pools (q) (0.31 ≤ λ ≤ λ0.45) and microbial growth rates (I max ) (−0.55 ≤ λ ≤ 0.67).Biomass accumulation is moderately sensitive to death rates (α A and α H ), the efficiency of heterotrophic growth (Y H ) and the allochthonous substrate input (v Sub ) (±0.15 ≤ λ ≤ ±0.41).Biomass accumulation has relatively low sensitivity (λ ≤ ±0.15) to variation in half-saturation constants (K L , K S , K N and K P ), parameters affecting only the dynamics of subglacial microbes (A 1 and H 1 ) (p Sub and k Sub ) and nitrogen fixers (A 3 and H 3 ) (n f , K N2 ), the bioavailability of labile substrate (J S1 ), exudate rates (ex A and ex H ), the input of allochthonous nutrients (v DIN and v DIP ) and the efficiency of autotrophic growth (Y A ). Variation in the half saturation for nitrogen (K N ) and phosphorus (K P ) causes little change to the accumulation of biomass (−0.02 ≤ λ ≤ 0.00 and −0.05 ≤ λ ≤ 0.03, respectively), but has a proportionally large effect on the accumulation of DIN, DIP and total nitrogen fixation (0.22 ≤ λ ≤ 0.95).Similarly, the reduction of efficiency and growth rates (n f ) for nitrogen fixers whilst fixing nitrogen (rather than assimilating DIN) has a relatively minor effect on the accumulation of biomass (λ ≤ ±0.09) but strongly affects DIN, DIP and total nitrogen fixed (λ = 0.96, λ = −0.96and λ = 0.60, respectively).
Microbial communities alter their metabolic state (through the Q 10 formulation) to persist during long periods of cold.At cold temperatures typical of glacier forefield environments, high Q 10 responses to temperature variation (Schipper et al., 2014) promote the survival of biomass under prolonged periods of harsh environmental conditions (soil temperatures < T ref ).
Heterotrophic production is critical in supporting the overall establishment of biomass and activity of the entire microbial community.Increasing the maximum growth rate of heterotrophs (I maxH ) leads to a substantial increase in both heterotrophic and autotrophic biomass (λ = 0.47 and λ = 0.67, respectively).However, strikingly, an increase in the maximum growth rate of autotrophs (I maxA ) has the net effect of lowering autotrophic and heterotrophic biomass (λ = −0.36 and λ = −0.55,respectively).Autotrophic growth is responsible for the build-up of substrate during the initial stages of soil development, thus supporting heterotrophic production.However, the dominance of autotrophic communities rapidly consumes available nutrients, at the expense of heterotrophs.If autotrophs are given a competitive advantage (i.e.increasing I maxA ), their rapid growth increases nutrient scarcity (λ DIN = −0.57,λ DIP = −1.58)and heterotrophic growth becomes nutrient limited.However, if heterotrophs are given a competitive advantage (i.e.increasing I maxH ), substrate is degraded more rapidly, liberating DIN (λ = 0.76) and DIP (λ = 0.38).This is an effective positive feedback effect, whereby the additional nutrients recycled from rapid heterotrophic degradation are able to support the growth of all microbial populations (including autotrophs (λ = 0.47)) in this oligotrophic and relatively nutrient poor environment.Similarly, increasing autotrophic death rates (α A ) reduces competition for nutrient resources whilst also increasing the availability of degradable organic matter (autotrophic necromass), thereby increasing heterotrophic biomass (λ = 0.24).
An increase in the heterotrophic growth efficiency (Y H ) leads to an overall more rapid accumulation of heterotrophic biomass (λ = 0.14) at the expense of autotrophs (λ = −0.41),since heterotrophic nutrient uptake is higher.A decrease in heterotrophic growth efficiency effectively means that for heterotrophs to grow by the same amount, they must degrade more organic matter per mole of carbon incorporated in biomass, thereby liberating nutrients and in turn supporting autotrophic growth.Previous modelling studies have consistently found microbial growth efficiency (Y ) to be the most sensitive parameter in determining overall biomass accumulation and the compartmentalisation of carbon between biotic and abiotic pools (Blagodatsky and Richter, 1998;Ingwersen et al., 2008;Toal et al., 2000).

Linearity
Linearity of parameter sensitivity is explored qualitatively over the range of parameter values by plotting the change in model output and visually fitting a linear or polynomial relationship to the trend (Fig. 9).Sensitivity and linearity evaluation have important implications in model optimisation.The highest degrees of freedom can be given to those parameters with fairly low sensitivity and a high degree of linearity (Q 10 , ex A , ex H , K L , K S , K N and K P ), since these parameters affect model output minimally and in a reasonably predictable way (illustrated in Fig. 6a).The parameters v Sub , p Sub , q and Y H also behave fairly linearly.Changes in model output respond non-linearly to changes in I maxA , I maxH , J S2 , d, α A , α H , v DIN and v DIP .Therefore, the sensitivity value λ is likely to change when looking at values away from the nominal parameter value, as illustrated for I maxA in Fig. 6b.For nonlinear parameters α A , α H , v DIN and v DIP , the maximum sensitivity is found in the region near the nominal value.This means that slight changes in parameter values will greatly affect the model output.Parameters I maxA , I maxH , J S2 , and d behave fairly linearly with comparatively little sensitivity around the nominal value; however their sensitivity increases with distance from the nominal value.This may give the impression of stability around the nominal values however a tipping point may be reached when parameters deviate too much from the nominal value.As such, non-linear parameters (identified in Fig. 9) must be given due caution in optimisation exercises since changes in these parameters may yield unexpected model behaviour.

Uncertainty
Uncertainty is evaluated over the entire plausible parameter range and results are presented in Fig. 10.The parameters that bear the highest uncertainty in the accumulation of autotrophic and heterotrophic biomass are I maxH (397 and 685 %, respectively) and J S2 (333 and 606 %, respectively).A high degree of uncertainty (≥ 60 %) also results from variation in α A , α H , I maxA , q, Y H , d, v Sub and v DIP .This is due to a combination of high sensitivity and large plausible ranges.
There is minimal uncertainty (≤ 30 %) resulting from variation in parameters ex A , ex H , p Sub , k Sub , K L , K N , K P , K N2 , J S1 , Y A and v DIN .Parameters that bear high uncertainty in the accumulation of biomass also tend to cause high degrees of uncertainty in other model outputs, most notably the accumulation of DIN and DIP, total nitrogen fixed and seasonality.
Measurements of bacterial growth are fundamental to most aspects of microbial ecology.Consequently, there are many estimates for I maxA and I maxH from literature and related modelling studies; however they span over an order of magnitude (0.24 to 4.80 d −1 ) (Mur et al., 1999;Van Liere and Walsby, 1982;Frey et al., 2010;Ingwersen et al., 2008;Knapp et al., 1983;Zelenev et al., 2000;Stapleton et al., 2005; Darrah, 1991;Blagodatsky and Richter, 1998;Vandewerf and Verstraete, 1987a;Foereid and Yearsley, 2004;Toal et al., 2000;Scott et al., 1995), greatly increasing the uncertainty associated with I maxA and I maxH (average uncertainty = 265 and 693 %, respectively).Maximum growth rates have been experimentally measured for soils from the Damma Glacier, with a value roughly in the middle of the plausible range (Frey et al., 2010); however there is inherent abstraction when incorporating laboratory measurements  into models because of the assumptions and simplifications in model design.

Non-linear & most sensitive at nominal parameter values
Microbial death rates (α A and α H ), however, are difficult and problematic to define experimentally (Toal et al., 2000), and there is a great deal of variation in how losses from microbial biomass are modelled.Death rates bear moderate sensitivity in model outputs (λ ≤ −0.13 and λ ≥ 0.24), and very high uncertainty (77 to 178 %), as a consequence of the large plausible parameter range, which spans almost 2 orders of magnitude (0.006 to 0.48 d −1 ).Furthermore, the transferability of microbial death rate constants is compromised by the different mathematical formulations used to describe death rates (e.g.constant fixed rate (German et al., 2012;Grant et al., 1993;Scott et al., 1995;Toal et al., 2000): logistic (Boudreau, 1999;Kravchenko et al., 2004) and variable, depending on model conditions (Knapp et al., 1983;Blagodatsky and Richter, 1998;Ingwersen et al., 2008;Zelenev et al., 2000;Lancelot et al., 2005)).
The allochthonous deposition of organic matter and nutrients is a source of considerable uncertainty.Published literature provides estimates of snow nutrient concentrations at the Damma Glacier (Brankatschk et al., 2011); however the fate of these nutrients once deposited on the soil surface is largely unknown.Accordingly, the plausible parameter ranges for v Sub , v DIN and v DIP are wide, and the resulting uncertainty is very high (75 to 2503 %).Furthermore, variation in parameters v DIN and v DIP results in highly non-linear behaviour, with comparatively large changes resulting from small changes in the parameter value.However, the effect of v DIN on the accumulation of autotrophic and heterotrophic biomass is minimal (5 and 9 %), due to the substantial inputs of nitrogen from nitrogen fixers throughout the development of the forefield.Analysis of the BACWAVE model (Zelenev et al., 2000) found high sensitivity of the spatial and temporal response of bacterial populations to changes in allochthonous carbon sources in soil.
The partitioning of organic matter compounds into a limited number of substrate pools (e.g.labile and refractory) is, although common practice, artificial, and a meaningful value cannot be determined experimentally or from literature (e.g.Arndt et al., 2013).Therefore, these parameters are ideal for initial tuning and optimisation exercises.Adjustment in q, J S1 and J S2 causes considerable variation in the accumulation of autotrophic and heterotrophic biomass (0.07 ≤ λ ≤ 0.64) and relatively large uncertainty (4 to 606 %).
Estimation of the inhibition of nitrogen fixation with DIN (K N2 ) is based on published literature and a relatively large range is considered in the uncertainty analysis (56.4 µg N g −1 ≤ K N2 ≤ 1128 µg N g −1 ) (LaRoche and Breitbarth, 2005;Rabouille et al., 2006;Holl and Montoya, 2005).The total nitrogen fixed and total accumulation of DIN is fairly sensitive to variation in K N2 (λ = 0.08 and λ = 0.11, respectively), however the resulting uncertainty in the overall accumulation of autotrophic and heterotrophic biomass is low (3 and 1 %, respectively).Therefore, although this parameter is not well defined, its importance is outweighed by other much more sensitive parameters such as the active fraction (d) and microbial death rates (α A and α H ).
Many of the fairly well-constrained parameters result in low uncertainty values.This can be explained by relatively tight parameter bounds explored in the uncertainty analysis, or relatively low sensitivity of model output to variation in this parameter, or a combination of both of these factors.For example, exudate production rates (ex A and ex H ) (uncertainty ≤ 2 %) have relatively tight parameter bounds (Allison, 2005), as well as low overall sensitivity (−0.01 ≤ λ ≤ 0.02).Half-saturation constants for carbon substrate (K S ) have been estimated experimentally (Vandewerf and Verstraete, 1987a;Vandewerf and Verstraete, 1987b;Blagodatsky et al., 1998;Anderson and Domsch, 1985) and fitted using models (Darrah, 1991;Ingwersen et al., 2008;Stapleton et al., 2005), with values between 50 and 1000 µg C g −1 .This variation can be attributed to differences in experimental technique and medium of substrate used, and differences in model structure and optimisation; however overall model sensitivity is low (λ = −0.07),as is reflected in other models (Ingwersen et al., 2008;Blagodatsky and Richter, 1998).

Case study 1: Damma Glacier
Model parameters are optimised to obtain the best possible fit to observational data from the Damma Glacier, Switzerland (Table 5 and Fig. 11).Total microbial biomass increases rapidly during the initial stages of the soil development to a peak of 61.6 µg C g −1 (year 12), followed by a decline over the following 10 years, after which biomass is fairly stable at roughly 37.0 µg C g −1 (years 30 to 75).The microbial com-munity evolves from an even split between all six pools of microbial biomass (16.7 %) to a community dominated by autotrophs (A 1−3 comprises 78.5 % of biomass in year 13) (Fig. 12).Nitrogen-fixing autotrophs (A 3 ) are the dominant functional group during the first 10 years of the simulation (up to 56.1 % of biomass), after which soil autotrophs (A 2 ) increase in relative abundance (up to 35.6 % of biomass) as DIN concentrations increase.Subglacial microbes (A 1 and H 1 ) are consistently outcompeted by soil microbes (A 2 and H 2 ) (Fig. 12).Primary production (A 1−3 ) accounts for between 68.7 and 88.9 % of total bacterial production, whereas heterotrophic production (H 1−3 ) accounts for the remaining 11.1 to 31.3 % (Fig. 13b).This trend is also reflected in independent field studies, whereby autotrophic production has been identified as a major source of carbon in young soils at the Damma Glacier (Zumsteg et al., 2013b;Esperschütz et al., 2011;Frey et al., 2013) and elsewhere including the Puca Glacier in Peru (Schmidt et al., 2008).Heterotrophic production is closely associated with the abundance and availability of carbon substrate.A high proportion of labile substrate (39.4 % in year 1) supports high rates of heterotrophic production and rapid accumulation of heterotrophic biomass.rapidly depleted (Fig. 13a) followed by a sharp decline in biomass (Fig. 11).Following the exhaustion of labile organic carbon substrate, heterotrophic production is sustained at lower steady rate (roughly 10.0 µg C g −1 yr −1 ) and predicted microbial biomass is within the natural variability of the observational data.Chemical analysis of substrate from the Damma Glacier forefield suggests that organic matter becomes increasingly refractory in the later stages of development due to continual re-working and cycling by microbial communities (Goransson et al., 2011), as reflected in the model.Soil respiration (net DIC efflux) follows a broadly similar pattern to total microbial production, and is relatively stable (roughly 312.0 µg C g −1 yr −1 ) in the later stages of soil development (years 30 to 75).Soil respiration rates in the Damma Glacier have been estimated to be in the range of 130.0 µg C g −1 yr −1 (Schulz et al., 2013;Guelland et al., 2013a), which is within the range predicted by the SHIM-MER model.
No parameter combination could reproduce the high substrate accumulation observed at the Damma Glacier (Fig. 11).Even under extremely high allochthonous substrate loading (v Sub = 3.0, equivalent of 195.5 µg C g −1 yr −1 ), carbon substrate accumulates to roughly half of the highest maximum substrate sampled in the BigLink project (31 363.1 µg C g −1 ) (Bernasconi et al., 2011).We attribute this to the extremely rapid onset of vegetation (Bernasconi et al., 2011).Duc et al. (2009) compare rhizosphere and bulk soils in the Damma Glacier, and find substantially higher total organic carbon concentrations in soils sampled in close proximity to plants.The SHIMMER model does not include a vegetation component and is thus not able to account for the effect of plants.
Field-based nutrient enrichment experiments show that microbial growth is limited by carbon and nitrogen (Gorans-son et al., 2011;Yoshitake et al., 2007).Furthermore, a high diversity of diazotrophs has been associated with soil nitrogen accumulation in initial soils at the Damma Glacier (Duc et al., 2009), the Puca Glacier (Peru) (Schmidt et al., 2008;Nemergut et al., 2007), Mendenhall Glacier (Alaska) (Sattin et al., 2009;Knelman et al., 2012) and Anvers Island (Antarctica) (Strauss et al., 2012).This is reflected in the model.DIN is the primary limiting nutrient for subglacial and soil species (A 1 , A 2 , H 1 and H 2 ) during the first 5 years of simulated soil development, after which there is sufficient DIN in the soil (6.22 µg N g −1 ) to support the net accumulation of microbial biomass in all pools.However, nitrogen fixers (A 3 and H 3 ) are able to alleviate DIN limitation and experience net growth immediately, contributing 85.7 % of all nitrogen assimilated in microbial biomass during the first 10 years (Fig. 13c).The main supply of phosphorus to natural terrestrial ecosystems is the underlying parent rock, especially in glaciated settings whereby relatively high erosion rates and crushed rock flour give rise to increased mineral dissolution rates.However, phosphate is a relatively immobile macronutrient due to sorption and interaction with other soil constituents, making it a potential growth-limiting nutrient in terrestrial ecosystems (Hinsinger, 2001).Initial phosphorus limitation is rapidly alleviated by the accumulation of DIP at an average rate of 2.0 µg P g −1 yr −1 (years 1 to 10).Predicted DIP closely resembles field data (Goransson et al., 2011).Isotopic analysis (Tamburini et al., 2012) and modelling suggests that biological activity is the main driver of phosphorus cycling in developing soils at the Damma Glacier.
The seasonal evolution of the Damma Glacier forefield is not well understood; however transplantation studies have indicated that microbial communities respond dynamically to changing environmental conditions (Zumsteg et al., 2013a) and soil bacteria photosynthesise, degrade organic material and fix nitrogen at varying rates over spring, summer and autumn (Lazzaro et al., 2012).Microbial abundance and production calculated over winter (January, February, March) and summer (July, August, September) varies seasonally (Fig. 11a).Snow cover attenuates PAR in the winter period causing large seasonal fluctuations in the biomass and overall production of photosynthetic organisms (A 2 and A 3 ).Whilst heterotrophic production is also higher in the summer (0.13 µg C g −1 d −1 ) than winter (0.05 µg C g −1 d −1 ), populations remain more stable.SHIMMER estimates that 69.6 to 74.5 % of total net CO 2 efflux occurs during the 4 month snow-free summer period between June and October, which agrees well with estimations from field studies at the Damma Glacier (62 to 70 %) (Guelland et al., 2013b).

Case study 2: Athabasca Glacier, Canada
The model is calibrated and validated against a second test data set from the Athabasca Glacier forefield in Canada.Parameters (Table 5) are varied sequentially to provide a best fit to observations of microbial biomass and carbon substrate (Fig. 14).Microbial biomass accumulates throughout the simulation from 2.6 µg C g −1 in year 1 to 11.5 µg C g −1 in year 40, and at roughly 0.02 µg C g −1 yr −1 in years 50 to 75.These results agree with the observed accumulation of microbial biomass in the validation data set.During the early stages of succession (< 15 years), nitrogen-fixing autotrophs account for up to 62.5 % of the total biomass.Heterotrophs (H  (11.3 µg N g −1 ); however later increase to a similar relative abundance to nitrogen-fixing heterotrophs when DIN stocks are more plentiful.The model predicts the observed accumulation of carbon substrate with a substantially lower allochthonous input of substrate (v Sub = 0.05) compared to the Damma Glacier (v Sub = 0.60).Autotrophic production accounts for > 87.5 % of total bacterial production.DIN remains low compared to the Damma Glacier.The majority of the DIN assimilation in young soils (< 15 years) is by nitrogen fixation (up to 96.6 %).The seasonal oscillations in microbial biomass and activity at the Athabasca Glacier forefield are considerably smaller than the Damma Glacier forefield, due to increased nutrient scarcity (inhibiting growth and slowing the biotic response to seasonal variability) and lower microbial biomass.A high level of bacterial production is sustained by a continuous pool of labile substrate.Our ability to put these model results into context with field data is limited since there the Athabasca Glacier forefield is considerably less intensively studied than the Damma Glacier.Nevertheless, the Athabasca Glacier site acts as a secondary test and validation of model behaviour against published obser-vational and experimental data suggesting an increase in microbial biomass and accumulation of organic carbon.

Model development, application and recommendations
The SHIMMER framework greatly improves our ability to quantify dominant biogeochemical processes and make scenario-based predictions of soil development.An additional result of developing the SHIMMER framework has been an appreciation of the types of data that are necessary to build a mechanistic and fully constrained numerical representation of these systems, and this is presented here as a recommendation to future field and laboratory efforts.

Data availability
Model development and confidence in model evaluations would be improved by higher temporal and spatial sampling resolution along the chronosequence and interdisciplinary (microbial and geochemical) data sets (concentrations and rates).11).Therefore, there is a need for test sites where the influence of vegetation is not so great.A full mass budget of inputs (e.g.aerial deposition) and outputs (e.g.leaching) of substrate and nutrients through the soil surface would greatly improve the predictive power of SHIM-MER, as demonstrated in sensitivity analysis (Sect.5.2.1).

Death and dormancy
Many biological processes including microbial growth and death are simplified in order to describe them mathematically.Cell death is not well constrained, with a lack of empirical measurements, and fundamental differences in how these rates are defined in models (German et al., 2012;Grant et al., 1993;Scott et al., 1995;Toal et al., 2000;Boudreau, 1999;Kravchenko et al., 2004;Knapp et al., 1983;Blagodatsky and Richter, 1998;Ingwersen et al., 2008;Zelenev et al., 2000;Lancelot et al., 2005).Additionally, in situ measurements of active and dormant microbes in field studies are rare, thus making validation difficult.Those that are provided in the literature (Toal et al., 2000;Blagodatsky and Richter, 1998;Vandewerf and Verstraete, 1987b) are mostly derived from models and are a snap shot in time.These rates strongly influence both the microbial populations' stability in size, and the available substrate and necromass.Therefore, more focussed empirical observations are needed to support and inform model processes, increase confidence in predictions and support future development.

Bacterial growth
Parameters defined as constants in mathematical models are often known to vary in time depending on prevailing environmental and biogeochemical conditions, for example sensitivity to temperature, soil moisture, oxygen availability, C : N ratio and quality of soil organic matter (Manzoni et al., 2012;Erhagen et al., 2015).Variable growth efficiencies and background maintenance, for example, would require additional parameters that currently cannot be defined with sufficient confidence, increasing uncertainty rather than improving the model.Laboratory incubations are useful to estimate parameter values experimentally (Blagodatsky et al., 1998).It is our intention that future laboratory analysis can at least in part inform parameters such that confidence in model output is increased.

Discrete pools vs. continuum representation
In a real soil system, the composition of the microbial community, along with the quality of substrate are continua, rather than discreet categories (as in SHIMMER).Classifying variables into categories (e.g.labile and refractory) ultimately provides a simplistic view of the system, but also provides tuning parameters, flexibility and a high degree of generality.Other modelling approaches such as individualbased modelling (IBM) would account for population heterogeneity, and the ability to link mechanisms to population dynamics at the individual level (Hellweger and Bucci, 2009).However, this is outside the scope of the current model version.

Deterministic vs. stochastic
The biogeochemical development of a de-glaciated forefield is extremely heterogeneous across a range of spatial scales, which affects the signal of biogeochemical changes in observations (Bernasconi et al., 2011;Bradley et al., 2014).Stochastic variability resulting from disturbances, heterogeneity in biotic and abiotic processes, lateral and vertical particle movement and diffusion, inter-annual variability, and changing environmental conditions affect the biogeochemical development of a chronosequence in such a way that a deterministic model cannot predict.Furthermore, the SHIM-MER model does not predict or account for the effect of plant biomass, which is highly abundant in the forefield of the Damma Glacier (Bernasconi et al., 2011).However, a first attempt to gain quantitative insight into glacier forefield dynamics favours a deterministic framework.

0-D vs. multi-D
Currently, nearly all soil system models assume that at the finest level of detail, soil is a well-mixed homogeneous particle with respect to its composition and dynamics, whereas microbial populations and metabolic rates are known to be heterogeneous across a number of different spatial scales and directions.A large amount of additional parameters and equations would need to be incorporated in order to include 1-D or 2-D processes in SHIMMER.It is currently unfeasible to incorporate this level of detail due to limitations in the observational data, and the resulting model would not be useful for the purpose it is intended for, that is to describe, predict and provide insight on the development of initial ecosystems such as those exposed due to glacier retreat.

Conclusions
Accurate quantitative prediction of the biogeochemical development of de-glaciated forefields is important to understand the primary succession of microbial communities and the formation of organic carbon in extreme oligotrophic environments.The forefield ecosystem reacts rapidly to climate change (Smittenberg et al., 2012), and the fine glacial flour, highly reactive sediments and rapid biological cycling of nutrients typical of forefields may be significant to global biogeochemical cycles in the context of future large-scale ice retreat.
www.geosci-model-dev.net/8/3441/2015/Geosci.Model Dev., 8, 3441-3470, 2015 Here we present SHIMMER, a novel modelling framework designed to predict microbial community development during the initial stages of ecosystem formation.The model accurately predicts the accumulation of microbial biomass and organic carbon during the initial stages of soil development from two glacier forefields (Bernasconi et al., 2011;Goransson et al., 2011;Insam and Haselwandter, 1989), and supports our general understanding of these ecosystems.Autotrophic production and nitrogen fixation are fundamental to the establishment of microbial communities and stable and labile pools of organic substrate and inorganic nutrients, a finding that is supported by field experiments at the Damma Glacier (Zumsteg et al., 2013b;Esperschütz et al., 2011;Frey et al., 2013;Duc et al., 2009) and elsewhere including the Puca Glacier in Peru (Schmidt et al., 2008;Nemergut et al., 2007), the Mendenhall Glacier in Alaska (Sattin et al., 2009;Knelman et al., 2012), andAnvers Island in Antarctica (Strauss et al., 2012).Soil respiration is comparable to field observations (Schulz et al., 2013;Guelland et al., 2013a).The seasonal evolution of glacier forefields is not well understood (Bradley et al., 2014); however modelling work is likely to provide insight into the dynamics of the "non-growing season".For example, modelled summer production accounts for roughly 70 % of total annual respiration, in line with field observations from the Damma Glacier (Guelland et al., 2013b).
The accumulation of microbial biomass is highly sensitive to variation in the Q 10 values, active biomass, the bioavailability of organic matter, bacterial growth efficiency, and rates of microbial growth and death.These parameters also bear high uncertainty due to a relatively large range of plausible values.Many of the well-constrained parameters (e.g.half-saturation constants and exudation rates) have low sensitivity and uncertainty and show mostly linear behaviour.One of the striking outcomes of the sensitivity study is the apparent strong dependence between the heterotrophic and autotrophic microbial communities.Heterotrophic production degrades organic matter and recycles nutrients, in turn supporting autotrophic and heterotrophic growth.Increasing heterotrophic growth rates therefore increases the accumulation in all pools of biomass.
SHIMMER is the first step towards an iterative and interdisciplinary framework (presented in Fig. 2), integrating fieldwork and laboratory experiments with model development, testing and application.The development of SHIM-MER (1.0) is informed by previous experimental and field campaigns (Bradley et al., 2014).It is expected that further quantitative analysis of forefield dynamics will guide and inform future studies that will provide new data and insights, which will inform further model development and so forth.SHIMMER thus contributes to more accurate quantitative predictions that enable a deeper understanding of the processes which control microbial communities, their role on global biogeochemical cycles and their response to climate variations in the future.

Figure 3 .
Figure 3.A conceptual model showing the components and transfers of SHIMMER.State variables are indicated with shading.

-
total C substrate (average over the entire simulation run) -DIN (average over the entire simulation run) -DIP (average over the entire simulation run) -total ON (average over the entire simulation run) -total OP (average over the entire simulation run) -total nitrogen fixed (cumulative) -seasonal variation in microbial biomass (final year of simulation).

Figure 5 .
Figure 5. Illustration of calculation of sensitivity (λ) where (a) the value of λ is representative of the sensitivity; (b) the value of λ is not representative of the sensitivity.In (b) the apparent sensitivity (λ) will be low due to model behaviour either side of the nominal parameter value having an opposite sign, even though the model may be truly sensitive to that parameter.

Figure 6 .
Figure 6.Simulated response of autotrophic (A 1−3 ) biomass (blue line) and heterotrophic (H 1−3 ) biomass (red line) to variation in (a) Q 10 and (b) I maxH across the entire range of plausible values, forced with meteorological data from the Damma Glacier over 75 years.The shaded segment shows the region in which model sensitivity (λ) is calculated.

Figure 7 .
Figure 7. Model output of the Damma Glacier forefield for a 75 year time series with nominal parameter values (Table4).

Figure 8 .
Figure 8. Sensitivity of model outputs (λ) to individual parameter variation.The model is forced with meteorological data from the Damma Glacier (Fig. 4) over 75 years.

Figure 9 .
Figure 9. Heat map showing parameter linearity and non-linearity over a range of model outputs.The model is forced with meteorological data from the Damma Glacier (Fig. 4) over 75 years.

Figure 10 .
Figure 10.Heat map showing uncertainty of model outputs (ø) arising from individual parameters.The model is forced with meteorological data from the Damma Glacier (Fig. 4) over 75 years.

Figure 11 .
Figure 11.Model output optimised to observational data from the Damma Glacier, Switzerland.

Figure 13 .
Figure 13.Evolution of (a) substrate dynamics, (b) bacterial production and (c) nitrogen (N) assimilation of model output from the first 20 years of soil formation at the Damma Glacier, Switzerland.

Figure 14 .
Figure 14.Model output optimised to observational data from the Athabasca Glacier, Canada.

Table 2 .
State variables and initial values.
composition and net autotrophy/heterotrophy. SHIMMER does not explicitly account for vegetation and thus cannot reproduce the high organic carbon accumulation in vegetated sites (Fig.