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

J. Bradley, A. Anesio, J. Singarayer, M. Heath, S. Arndt

► To cite this version:

J. Bradley, A. Anesio, J. Singarayer, M. Heath, S. Arndt. SHIMMER (1.0): a novel mathematical model for microbial and biogeochemical dynamics in glacier forefield ecosystems. Geoscientific Model Development, 2015, 8 (10), pp.3441-3470. 10.5194/gmd-8-3441-2015. hal-04450531

HAL Id: hal-04450531

https://hal.science/hal-04450531

Submitted on 6 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

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

J. A. Bradley<fnref n="1" />,<fnref n="2" />, A. M. Anesio<fnref n="1" />, J. S. Singarayer<fnref n="3" />, M. R. Heath<fnref n="4" />, and S. Arndt<fnref n="2" />

*Correspondence to:* J. A. Bradley (j.bradley@bristol.ac.uk)

Received: 9 July 2015 – Published in Geosci. Model Dev. Discuss.: 6 August 2015

Revised: 15 October 2015 – Accepted: 17 October 2015 – Published: 28 October 2015

**Abstract.** SHIMMER (Soil biogeoChemIcal Model for Mi-crobial Ecosystem Response) is a new numerical modelling framework designed to simulate microbial dynamics and bio-geochemical cycling during initial ecosystem development in glacier forefield soils. However, it is also transferable to other extreme ecosystem types (such as desert soils or the sur-face of glaciers). The rationale for model development arises from decades of empirical observations in glacier forefields, and enables a quantitative and process focussed approach. Here, we provide a detailed description of SHIMMER, test its performance in two case study forefields: 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 accu-mulation of microbial biomass is highly dependent on vari-ation in microbial growth and death rate constants, Q₁₀ val-ues, the active fraction of microbial biomass and the reactiv-ity of organic matter. The model correctly predicts the rapid accumulation of microbial biomass observed during the ini-tial stages of succession in the forefields of both the case study systems. Primary production is responsible for the ini-tial build-up of labile substrate that subsequently supports heterotrophic growth. However, allochthonous contributions of organic matter, and nitrogen fixation, are important in sus-taining this productivity. The development and application of SHIMMER also highlights aspects of these systems that re-quire further empirical research: quantifying nutrient budgets and biogeochemical rates, exploring seasonality and micro-bial growth and cell death. This will lead to increased under-

standing of how glacier forefields contribute to global bio-geochemical cycling and climate under future ice retreat.

## 1 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<fnref n="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 thou-sands of years, are typically highly oligotrophic, with low nutrient budgets and are subject to harsh and rapidly fluctuat-ing environmental conditions. They potentially play a signif-icant yet largely unexplored role in large-scale biogeochemi-cal 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 phos-phorus cycle (Filippelli, 2002; Follmi et al., 2009) and the productivity of downstream and coastal ecosystems (Ane-sio et al., 2009; Mindl et al., 2007; Fountain et al., 2008; Anderson et al., 2000). Furthermore, the initial stages of

<footnote n="1">¹Bristol Glaciology Centre, School of Geographical Sciences, University of Bristol, Bristol, BS8 1SS, UK</footnote>

<footnote n="2">²BRIDGE, School of Geographical Sciences, University of Bristol, Bristol, BS8 1SS, UK</footnote>

<footnote n="3">³Department of Meteorology, University of Reading, Reading, RG6 6BB, UK</footnote>

<footnote n="4">⁴Department of Mathematics and Statistics, University of Strathclyde, Glasgow, G1 1XH, UK</footnote>
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

J. A. Bradley<fnref n="1" />,<fnref n="2" />, A. M. Anesio<fnref n="1" />, J. S. Singarayer<fnref n="3" />, M. R. Heath<fnref n="4" />, and S. Arndt<fnref n="2" />

*Correspondence to:* J. A. Bradley (j.bradley@bristol.ac.uk)

Received: 9 July 2015 – Published in Geosci. Model Dev. Discuss.: 6 August 2015

Revised: 15 October 2015 – Accepted: 17 October 2015 – Published: 28 October 2015

**Abstract.** SHIMMER (Soil biogeoChemIcal Model for Mi-crobial Ecosystem Response) is a new numerical modelling framework designed to simulate microbial dynamics and bio-geochemical cycling during initial ecosystem development in glacier forefield soils. However, it is also transferable to other extreme ecosystem types (such as desert soils or the sur-face of glaciers). The rationale for model development arises from decades of empirical observations in glacier forefields, and enables a quantitative and process focussed approach. Here, we provide a detailed description of SHIMMER, test its performance in two case study forefields: 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 accu-mulation of microbial biomass is highly dependent on vari-ation in microbial growth and death rate constants, Q₁₀ val-ues, the active fraction of microbial biomass and the reactiv-ity of organic matter. The model correctly predicts the rapid accumulation of microbial biomass observed during the ini-tial stages of succession in the forefields of both the case study systems. Primary production is responsible for the ini-tial build-up of labile substrate that subsequently supports heterotrophic growth. However, allochthonous contributions of organic matter, and nitrogen fixation, are important in sus-taining this productivity. The development and application of SHIMMER also highlights aspects of these systems that re-quire further empirical research: quantifying nutrient budgets and biogeochemical rates, exploring seasonality and micro-bial growth and cell death. This will lead to increased under-

standing of how glacier forefields contribute to global bio-geochemical cycling and climate under future ice retreat.

## 1 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<fnref n="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 thou-sands of years, are typically highly oligotrophic, with low nutrient budgets and are subject to harsh and rapidly fluctuat-ing environmental conditions. They potentially play a signif-icant yet largely unexplored role in large-scale biogeochemi-cal 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 phos-phorus cycle (Filippelli, 2002; Follmi et al., 2009) and the productivity of downstream and coastal ecosystems (Ane-sio et al., 2009; Mindl et al., 2007; Fountain et al., 2008; Anderson et al., 2000). Furthermore, the initial stages of

<footnote n="1">¹Bristol Glaciology Centre, School of Geographical Sciences, University of Bristol, Bristol, BS8 1SS, UK</footnote>

<footnote n="2">²BRIDGE, School of Geographical Sciences, University of Bristol, Bristol, BS8 1SS, UK</footnote>

<footnote n="3">³Department of Meteorology, University of Reading, Reading, RG6 6BB, UK</footnote>

<footnote n="4">⁴Department of Mathematics and Statistics, University of Strathclyde, Glasgow, G1 1XH, UK</footnote>
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.0 mg C g⁻¹, 0.1–2.0 mg N g⁻¹, and 2–8 µg P g⁻¹ (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: autotchthonous 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₂, 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, 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-
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.0 mg C g⁻¹, 0.1–2.0 mg N g⁻¹, and 2–8 µg P g⁻¹ (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₂, 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, 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><img src="image_3.png" /><figcaption>**Figure 1.** A conceptual overview of the main transformations and fluxes of nutrients and organic matter in and along a typical de-glaciated forefield.</figcaption></figure>

<table><caption>**Table 1.** Outline of quantitative and predictive modelling strategy.</caption><thead><tr><th>Quantitative modelling</th><th>Details</th></tr></thead><tbody><tr><td>Do initial microbial communities rely on autochthonous or allochthonous carbon?<br>Assess the role of nitrogen as a key limiting nutrient.</td><td>Explore scenarios of carbon loading to determine the reliance on allochthonous carbon vs. autochthonous production.<br>Quantify nitrogen budget to assess relative importance of nitrogen fixation vs. DIN assimilation in recently exposed and old soils.</td></tr><tr><td>Assess how microbial diversity affects soil development.</td><td>Test various community assemblages and characterise the soil environment they create.</td></tr><tr><th>Predictive modelling</th><th>Details</th></tr><tr><td>Identifying tipping points within the system.<br>Assess importance of seasonality.</td><td>Do population dynamics, external nutrient loading, seasonal changes, disturbance events and climate change trigger significant changes in soil development?<br>Explore the dynamics of the “non-growing season”, where there is little observational data due to inaccessibility and difficult field work conditions.</td></tr><tr><td>Assess importance of disturbance events and community reset.</td><td>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)</td></tr><tr><td>Assess sensitivity to future climate change</td><td>Explore the effect of climate and anthropogenic impact using scenario-based predictions.</td></tr><tr><td>Draw attention to gaps in our understanding and areas of future research.</td><td>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.</td></tr></tbody></table>

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. SHIMMER 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 SHIMMER 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,
<figure><img src="image_4.png" /><figcaption>**Figure 1.** A conceptual overview of the main transformations and fluxes of nutrients and organic matter in and along a typical de-glaciated forefield.</figcaption></figure>

<table><caption>**Table 1.** Outline of quantitative and predictive modelling strategy.</caption><thead><tr><th>Quantitative modelling</th><th>Details</th></tr></thead><tbody><tr><td>Do initial microbial communities rely on autochthonous or allochthonous carbon?<br>Assess the role of nitrogen as a key limiting nutrient.</td><td>Explore scenarios of carbon loading to determine the reliance on allochthonous carbon vs. autochthonous production.<br>Quantify nitrogen budget to assess relative importance of nitrogen fixation vs. DIN assimilation in recently exposed and old soils.</td></tr><tr><td>Assess how microbial diversity affects soil development.</td><td>Test various community assemblages and characterise the soil environment they create.</td></tr><tr><th>Predictive modelling</th><th>Details</th></tr><tr><td>Identifying tipping points within the system.<br>Assess importance of seasonality.</td><td>Do population dynamics, external nutrient loading, seasonal changes, disturbance events and climate change trigger significant changes in soil development?<br>Explore the dynamics of the “non-growing season”, where there is little observational data due to inaccessibility and difficult field work conditions.</td></tr><tr><td>Assess importance of disturbance events and community reset.</td><td>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)</td></tr><tr><td>Assess sensitivity to future climate change</td><td>Explore the effect of climate and anthropogenic impact using scenario-based predictions.</td></tr><tr><td>Draw attention to gaps in our understanding and areas of future research.</td><td>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.</td></tr></tbody></table>

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. SHIMMER 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 SHIMMER 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,
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Iterative, fully integrated interdisciplinary framework of SHIMMER, aimed at continuing model development as new data sets emerge and knowledge of the field evolves.</figcaption></figure>

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. SHIMMER 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.

## 2 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.

### 2.1 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. The model can be easily expanded into higher dimensions by including
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Iterative, fully integrated interdisciplinary framework of SHIMMER, aimed at continuing model development as new data sets emerge and knowledge of the field evolves.</figcaption></figure>

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. SHIMMER 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.

## 2 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.

### 2.1 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. The model can be easily expanded into higher dimensions by including
<figure><img src="image_5.png" /><figcaption>**Figure 3.** A conceptual model showing the components and transfers of SHIMMER. State variables are indicated with shading.</figcaption></figure>

<table><caption>**Table 2.** State variables and initial values.</caption><thead><tr><th>State variable</th><th>Units</th><th>Description</th><th>Damma Glacier, Switzerland</th><th>Reference</th><th>Athabasca Glacier, Canada</th><th>Reference</th></tr></thead><tbody><tr><td>A<sub>1</sub></td><td>µg C g<sup>-1</sup></td><td>Subglacial chemolithoautotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>A<sub>2</sub></td><td>µg C g<sup>-1</sup></td><td>Generic soil autotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>A<sub>3</sub></td><td>µg C g<sup>-1</sup></td><td>Nitrogen-fixing soil autotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>H<sub>1</sub></td><td>µg C g<sup>-1</sup></td><td>Subglacial heterotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>H<sub>2</sub></td><td>µg C g<sup>-1</sup></td><td>Generic soil heterotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>H<sub>3</sub></td><td>µg C g<sup>-1</sup></td><td>Nitrogen-fixing soil heterotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>S<sub>1</sub></td><td>µg C g<sup>-1</sup></td><td>Labile carbon substrate</td><td>278.520</td><td>Bernasconi et al. (2011)</td><td>117.720</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>S<sub>2</sub></td><td>µg C g<sup>-1</sup></td><td>Refractory carbon substrate</td><td>417.780</td><td>Bernasconi et al. (2011)</td><td>176.580</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>DIN</td><td>µg N g<sup>-1</sup></td><td>Dissolved inorganic nitrogen</td><td>0.160</td><td>Brankatschk et al. (2011)</td><td>0.070</td><td>Stoichiometric</td></tr><tr><td>DIP</td><td>µg P g<sup>-1</sup></td><td>Dissolved inorganic phosphorus</td><td>0.500</td><td>Bernasconi et al. (2011)</td><td>0.210</td><td>Stoichiometric</td></tr><tr><td>ON<sub>1</sub></td><td>µg N g<sup>-1</sup></td><td>Labile nitrogen substrate</td><td>39.440</td><td>Bernasconi et al. (2011)</td><td>16.600</td><td>Stoichiometric</td></tr><tr><td>ON<sub>2</sub></td><td>µg N g<sup>-1</sup></td><td>Refractory nitrogen substrate</td><td>59.160</td><td>Bernasconi et al. (2011)</td><td>24.900</td><td>Stoichiometric</td></tr><tr><td>OP<sub>1</sub></td><td>µg P g<sup>-1</sup></td><td>Labile phosphorus substrate</td><td>23.120</td><td>Stoichiometric</td><td>9.770</td><td>Stoichiometric</td></tr><tr><td>OP<sub>2</sub></td><td>µg P g<sup>-1</sup></td><td>Refractory phosphorus substrate</td><td>34.680</td><td>Stoichiometric</td><td>14.660</td><td>Stoichiometric</td></tr></tbody></table>

transport terms. The model represents the top centimetre of the soil surface as a homogeneous mix, and light, temperature, nutrients, organic compounds and microbial biomass are thus assumed to be evenly distributed. External nutrient inputs and leaching are additions to and subtractions from an external environment.

## 2.2 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).

### 2.2.1 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}$, $H_1-3$) (Table 2), to account for the different metabolic needs and physiological pathways of microbial organisms commonly found in glacier forefields, but remain manageable in terms of validating the model behaviour with existing data sets. $A_1$ and $H_1$ represent microbes derived residually from the subglacial environment, and will be referred to as “subglacial microbes”. $A_1$ represents chemolithoautotrophic microbes such as the genus *Thiobacillus*, and $H_1$ represents heterotrophic subglacial microbes such as the family Comamonadaceae. $A_2$ and $H_2$ represent microbes commonly found in glacier forefield soils with no “special” characteristics, such as *Leptolyngbya* and *Sphingopyxis*, respectively, and will be
<figure><img src="image_6.png" /><figcaption>**Figure 3.** A conceptual model showing the components and transfers of SHIMMER. State variables are indicated with shading.</figcaption></figure>

<table><caption>**Table 2.** State variables and initial values.</caption><thead><tr><th>State variable</th><th>Units</th><th>Description</th><th>Damma Glacier, Switzerland</th><th>Reference</th><th>Athabasca Glacier, Canada</th><th>Reference</th></tr></thead><tbody><tr><td>A<sub>1</sub></td><td>µg C g<sup>-1</sup></td><td>Subglacial chemolithoautotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>A<sub>2</sub></td><td>µg C g<sup>-1</sup></td><td>Generic soil autotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>A<sub>3</sub></td><td>µg C g<sup>-1</sup></td><td>Nitrogen-fixing soil autotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>H<sub>1</sub></td><td>µg C g<sup>-1</sup></td><td>Subglacial heterotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>H<sub>2</sub></td><td>µg C g<sup>-1</sup></td><td>Generic soil heterotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>H<sub>3</sub></td><td>µg C g<sup>-1</sup></td><td>Nitrogen-fixing soil heterotrophs</td><td>0.617</td><td>Bernasconi et al. (2011)</td><td>0.950</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>S<sub>1</sub></td><td>µg C g<sup>-1</sup></td><td>Labile carbon substrate</td><td>278.520</td><td>Bernasconi et al. (2011)</td><td>117.720</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>S<sub>2</sub></td><td>µg C g<sup>-1</sup></td><td>Refractory carbon substrate</td><td>417.780</td><td>Bernasconi et al. (2011)</td><td>176.580</td><td>Insam and Haselwandter (1989)</td></tr><tr><td>DIN</td><td>µg N g<sup>-1</sup></td><td>Dissolved inorganic nitrogen</td><td>0.160</td><td>Brankatschk et al. (2011)</td><td>0.070</td><td>Stoichiometric</td></tr><tr><td>DIP</td><td>µg P g<sup>-1</sup></td><td>Dissolved inorganic phosphorus</td><td>0.500</td><td>Bernasconi et al. (2011)</td><td>0.210</td><td>Stoichiometric</td></tr><tr><td>ON<sub>1</sub></td><td>µg N g<sup>-1</sup></td><td>Labile nitrogen substrate</td><td>39.440</td><td>Bernasconi et al. (2011)</td><td>16.600</td><td>Stoichiometric</td></tr><tr><td>ON<sub>2</sub></td><td>µg N g<sup>-1</sup></td><td>Refractory nitrogen substrate</td><td>59.160</td><td>Bernasconi et al. (2011)</td><td>24.900</td><td>Stoichiometric</td></tr><tr><td>OP<sub>1</sub></td><td>µg P g<sup>-1</sup></td><td>Labile phosphorus substrate</td><td>23.120</td><td>Stoichiometric</td><td>9.770</td><td>Stoichiometric</td></tr><tr><td>OP<sub>2</sub></td><td>µg P g<sup>-1</sup></td><td>Refractory phosphorus substrate</td><td>34.680</td><td>Stoichiometric</td><td>14.660</td><td>Stoichiometric</td></tr></tbody></table>

transport terms. The model represents the top centimetre of the soil surface as a homogeneous mix, and light, temperature, nutrients, organic compounds and microbial biomass are thus assumed to be evenly distributed. External nutrient inputs and leaching are additions to and subtractions from an external environment.

## 2.2 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).

### 2.2.1 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}$, $H_1-3$) (Table 2), to account for the different metabolic needs and physiological pathways of microbial organisms commonly found in glacier forefields, but remain manageable in terms of validating the model behaviour with existing data sets. $A_1$ and $H_1$ represent microbes derived residually from the subglacial environment, and will be referred to as “subglacial microbes”. $A_1$ represents chemolithoautotrophic microbes such as the genus *Thiobacillus*, and $H_1$ represents heterotrophic subglacial microbes such as the family Comamonadaceae. $A_2$ and $H_2$ represent microbes commonly found in glacier forefield soils with no “special” characteristics, such as *Leptolyngbya* and *Sphingopyxis*, respectively, and will be
<table><caption>**Table 3.** Derived variables.</caption><thead><tr><th>Derived variable</th><th>Description</th><th>Formulation</th></tr></thead><tbody><tr><td>cum_A<sub>1</sub></td><td>Increase in biomass of A<sub>1</sub></td><td>U<sub>A1</sub></td></tr><tr><td>cum_A<sub>2</sub></td><td>Increase in biomass of A<sub>2</sub></td><td>U<sub>A2</sub></td></tr><tr><td>cum_A<sub>3</sub></td><td>Increase in biomass of A<sub>3</sub></td><td>U<sub>A3</sub></td></tr><tr><td>cum_H<sub>1</sub></td><td>Increase in biomass of H<sub>1</sub></td><td>U<sub>H1</sub></td></tr><tr><td>cum_H<sub>2</sub></td><td>Increase in biomass of H<sub>2</sub></td><td>U<sub>H2</sub></td></tr><tr><td>cum_H<sub>3</sub></td><td>Increase in biomass of H<sub>3</sub></td><td>U<sub>H3</sub></td></tr><tr><td>cum_A</td><td>Increase in autotrophic biomass</td><td>ΣU<sub>A1-3</sub></td></tr><tr><td>cum_H</td><td>Increase in heterotrophic biomass</td><td>ΣU<sub>H1-3</sub></td></tr><tr><td>cum雄BP</td><td>Bacterial production</td><td>ΣU<sub>A1-3</sub> + ΣU<sub>H1-3</sub></td></tr><tr><td>cum_DIC_A</td><td>DIC produced by autotrophs</td><td>(1 - Y<sub>A</sub>) ⋅ (1/γ<sub>A</sub>) ⋅ U<sub>A1</sub> + (1 - Y<sub>A</sub>) ⋅ (1/γ<sub>A</sub>) ⋅ U<sub>A2</sub> + (1 - (Y<sub>A</sub> ⋅ n<sub>f</sub>)) ⋅ (1/γ<sub>A</sub> ⋅ n<sub>f</sub>) ⋅ U<sub>A3</sub></td></tr><tr><td>cum_DIC_H</td><td>DIC produced by heterotrophs</td><td>(1 - Y<sub>H</sub>) ⋅ (1/γ<sub>H</sub>) ⋅ U<sub>H1</sub> + (1 - Y<sub>H</sub>) ⋅ (1/γ<sub>H</sub>) ⋅ U<sub>H2</sub> + (1 - (Y<sub>H</sub> ⋅ n<sub>f</sub>)) ⋅ (1/γ<sub>H</sub> ⋅ n<sub>f</sub>) ⋅ U<sub>H3</sub></td></tr><tr><td>cum_DIC</td><td>Total DIC produced</td><td>cum_DIC_A + cum_DIC_H</td></tr><tr><td>cum_DIN</td><td>DIN consumed by all microbes</td><td>DIN<sub>Consumed</sub></td></tr><tr><td>cum_n<sub>f</sub></td><td>DIN fixed by A<sub>3</sub> and H<sub>3</sub></td><td>NC ⋅ ((1 - β) ⋅ U<sub>A3</sub> + (1 - β) ⋅ U<sub>H3</sub>)</td></tr></tbody></table>

referred to as “soil microbes”. A<sub>3</sub> and H<sub>3</sub> are “nitrogen-fixing microbes” or “nitrogen fixers”, able to fix atmospheric N<sub>2</sub> gas as a source of nitrogen when dissolved inorganic nitrogen (DIN) stocks become limiting, giving them a competitive advantage. Examples of A<sub>3</sub> and H<sub>3</sub> are *Nostoc* and *Rhizobiales*, respectively.

The overall rate of change of autotrophic and heterotrophic communities is described by

$$ \frac{dA_i}{dt} = U_{A_i} - G_{A_i} - X_{A_i}, \qquad (1) $$

$$ \frac{dH_i}{dt} = U_{H_i} - G_{H_i} - X_{H_i}. \qquad (2) $$

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<sub>A1</sub>), photo-autotrophic (U<sub>A2</sub> and U<sub>A3</sub>) and heterotrophic (U<sub>H1-3</sub>) growth. Subglacial chemolithoautotrophs (A<sub>1</sub>) acquire their energy through chemical synthesis of mineral compounds. Autotrophic soil species (A<sub>2</sub> and A<sub>3</sub>) 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

$$ U_{A1} = A_1 \cdot T_f \cdot d \cdot p_{\text{sub}} \cdot I_{\text{maxA}} \cdot \left( \frac{\text{DIN}}{\text{DIN} + (K_n \cdot k_{\text{Sub}})} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + (K_p \cdot k_{\text{Sub}})} \right), \qquad (3) $$

$$ U_{A2} = A_2 \cdot T_f \cdot d \cdot I_{\text{maxA}} \cdot \left( \frac{\text{PAR}}{\text{PAR} + K_L} \right) \cdot \left( \frac{\text{DIN}}{\text{DIN} + K_n} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + K_p} \right), \qquad (4) $$

where autotrophic growth is calculated according to a maximum growth rate (I<sub>max</sub>) (at temperature (T) = T<sub>ref</sub>) 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<sub>f</sub>) with a Q<sub>10</sub> formulation (Table 4) effectively slowing down or speeding up all life processes (Soetaert and Herman, 2009; Yoshitake et al., 2010; Schipper et al., 2014). T represents the temperature of the soil. The reference temperature (T<sub>ref</sub>) is the temperature at which rates (I<sub>max</sub> and α) are described. The response of growth rates to temperature allows microorganisms to contend with harsh environmental conditions and promotes their overall longevity in a transient natural setting. PAR sensitivity of photo-autotrophs (A<sub>2</sub> and A<sub>3</sub>) is described by Monod kinetics, and a half-saturation constant K<sub>L</sub> is fitted to adjust the growth rate with PAR. Phototrophic activity is limited at low irradiance, and inhibition at a very high irradiances is not represented in the SHIMMER framework for sake of maintaining a manageable level of complexity. Nutrient (DIN and dissolved inorganic phosphorus (DIP)) limitation is also expressed as Monod kinetics, with half-saturation (K<sub>N</sub> and K<sub>P</sub>) expressions. Subglacial
<table><caption>**Table 3.** Derived variables.</caption><thead><tr><th>Derived variable</th><th>Description</th><th>Formulation</th></tr></thead><tbody><tr><td>cum_A<sub>1</sub></td><td>Increase in biomass of A<sub>1</sub></td><td>U<sub>A1</sub></td></tr><tr><td>cum_A<sub>2</sub></td><td>Increase in biomass of A<sub>2</sub></td><td>U<sub>A2</sub></td></tr><tr><td>cum_A<sub>3</sub></td><td>Increase in biomass of A<sub>3</sub></td><td>U<sub>A3</sub></td></tr><tr><td>cum_H<sub>1</sub></td><td>Increase in biomass of H<sub>1</sub></td><td>U<sub>H1</sub></td></tr><tr><td>cum_H<sub>2</sub></td><td>Increase in biomass of H<sub>2</sub></td><td>U<sub>H2</sub></td></tr><tr><td>cum_H<sub>3</sub></td><td>Increase in biomass of H<sub>3</sub></td><td>U<sub>H3</sub></td></tr><tr><td>cum_A</td><td>Increase in autotrophic biomass</td><td>ΣU<sub>A1-3</sub></td></tr><tr><td>cum_H</td><td>Increase in heterotrophic biomass</td><td>ΣU<sub>H1-3</sub></td></tr><tr><td>cum_BP</td><td>Bacterial production</td><td>ΣU<sub>A1-3</sub> + ΣU<sub>H1-3</sub></td></tr><tr><td>cum_DIC_A</td><td>DIC produced by autotrophs</td><td>(1 - Y<sub>A</sub>) · <sup>1</sup>&frasl;<sub>Y<sub>A</sub></sub> · U<sub>A1</sub> + (1 - Y<sub>A</sub>) · <sup>1</sup>&frasl;<sub>Y<sub>A</sub></sub> · U<sub>A2</sub> + (1 - (Y<sub>A</sub> · n<sub>f</sub>)) · <sup>1</sup>&frasl;<sub>Y<sub>A</sub> · n<sub>f</sub></sub> · U<sub>A3</sub></td></tr><tr><td>cum_DIC_H</td><td>DIC produced by heterotrophs</td><td>(1 - Y<sub>H</sub>) · <sup>1</sup>&frasl;<sub>Y<sub>H</sub></sub> · U<sub>H1</sub> + (1 - Y<sub>H</sub>) · <sup>1</sup>&frasl;<sub>Y<sub>H</sub></sub> · U<sub>H2</sub> + (1 - (Y<sub>H</sub> · n<sub>f</sub>)) · <sup>1</sup>&frasl;<sub>Y<sub>H</sub> · n<sub>f</sub></sub> · U<sub>H3</sub></td></tr><tr><td>cum_DIC</td><td>Total DIC produced</td><td>cum_DIC_A + cum_DIC_H</td></tr><tr><td>cum_DIN</td><td>DIN consumed by all microbes</td><td>DIN<sub>Consumed</sub></td></tr><tr><td>cum_n<sub>f</sub></td><td>DIN fixed by A<sub>3</sub> and H<sub>3</sub></td><td>NC · ((1 - β) · U<sub>A3</sub> + (1 - β) · U<sub>H3</sub>)</td></tr></tbody></table>

referred to as “soil microbes”. A<sub>3</sub> and H<sub>3</sub> are “nitrogen-fixing microbes” or “nitrogen fixers”, able to fix atmospheric N<sub>2</sub> gas as a source of nitrogen when dissolved inorganic nitrogen (DIN) stocks become limiting, giving them a competitive advantage. Examples of A<sub>3</sub> and H<sub>3</sub> are *Nostoc* and *Rhizobiales*, respectively.

The overall rate of change of autotrophic and heterotrophic communities is described by

$$ \frac{dA_i}{dt} = U_{A_i} - G_{A_i} - X_{A_i}, \qquad (1) $$

$$ \frac{dH_i}{dt} = U_{H_i} - G_{H_i} - X_{H_i}. \qquad (2) $$

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<sub>A1</sub>), photo-autotrophic (U<sub>A2</sub> and U<sub>A3</sub>) and heterotrophic (U<sub>H1-3</sub>) growth. Subglacial chemolithoautotrophs (A<sub>1</sub>) acquire their energy through chemical synthesis of mineral compounds. Autotrophic soil species (A<sub>2</sub> and A<sub>3</sub>) 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

$$ U_{A1} = A_1 \cdot T_f \cdot d \cdot p_{\text{sub}} \cdot I_{\text{maxA}} \cdot \left( \frac{\text{DIN}}{\text{DIN} + (K_n \cdot k_{\text{Sub}})} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + (K_p \cdot k_{\text{Sub}})} \right), \qquad (3) $$

$$ U_{A2} = A_2 \cdot T_f \cdot d \cdot I_{\text{maxA}} \cdot \left( \frac{\text{PAR}}{\text{PAR} + K_L} \right) \cdot \left( \frac{\text{DIN}}{\text{DIN} + K_n} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + K_p} \right), \qquad (4) $$

where autotrophic growth is calculated according to a maximum growth rate (I<sub>max</sub>) (at temperature (T) = T<sub>ref</sub>) 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<sub>f</sub>) with a Q<sub>10</sub> formulation (Table 4) effectively slowing down or speeding up all life processes (Soetaert and Herman, 2009; Yoshitake et al., 2010; Schipper et al., 2014). T represents the temperature of the soil. The reference temperature (T<sub>ref</sub>) is the temperature at which rates (I<sub>max</sub> and α) are described. The response of growth rates to temperature allows microorganisms to contend with harsh environmental conditions and promotes their overall longevity in a transient natural setting. PAR sensitivity of photo-autotrophs (A<sub>2</sub> and A<sub>3</sub>) is described by Monod kinetics, and a half-saturation constant K<sub>L</sub> is fitted to adjust the growth rate with PAR. Phototrophic activity is limited at low irradiance, and inhibition at a very high irradiances is not represented in the SHIMMER framework for sake of maintaining a manageable level of complexity. Nutrient (DIN and dissolved inorganic phosphorus (DIP)) limitation is also expressed as Monod kinetics, with half-saturation (K<sub>N</sub> and K<sub>P</sub>) expressions. Subglacial
<table><caption>**Table 4.** SHIMMER model formulation.</caption><tbody><tr><td colspan="2"><strong>Fundamental balance equations</strong></td></tr><tr><td>Rate of change of autotrophic biomass (A<sub>1-3</sub>)</td><td>$\frac{dA_i}{dt} = U_{Ai} - G_{Ai} - X_{Ai}$</td></tr><tr><td>Rate of change of heterotrophic biomass (H<sub>1-3</sub>)</td><td>$\frac{dH_i}{dt} = U_{Hi} - G_{Hi} - X_{Hi}$</td></tr><tr><td>Rate of change of labile carbon substrate (S<sub>1</sub>)</td><td>$\frac{dS_1}{dt} = v_{Sub} \cdot I_{S1} + q \cdot G_{Ai} + q \cdot G_{Hi} + X_A + X_H - U_{S1Hi} - W_{S1}$</td></tr><tr><td>Rate of change of refractory carbon substrate (S<sub>2</sub>)</td><td>$\frac{dS_2}{dt} = v_{CS2} \cdot I_{S2} + (1-q) \cdot G_{Ai} + (1-q) \cdot G_{Hi} - U_{S2Hi} - W_{S2}$</td></tr><tr><td>Rate of change of dissolved inorganic nitrogen (DIN)</td><td>$\frac{dDIN}{dt} = v_{DIN} \cdot I_{DIN} - DIN_{Consumed} + DIN_{Released} - W_{DIN}$</td></tr><tr><td>Rate of change of dissolved inorganic phosphorus (DIP)</td><td>$\frac{dDIP}{dt} = v_{DIP} \cdot I_{DIP} - DIP_{Consumed} + DIP_{Released} - W_{DIP}$</td></tr><tr><td>Rate of change of labile organic nitrogen (ON<sub>1</sub>)</td><td>$\frac{dON_1}{dt} = v_{Sub} \cdot I_{ON1} - ON_1Accumulation - ON_1Degraded - W_{ON1}$</td></tr><tr><td>Rate of change of refractory organic nitrogen (ON<sub>2</sub>)</td><td>$\frac{dON_2}{dt} = v_{Sub} \cdot I_{OP2} - OP_2Accumulation - OP_2Degraded - W_{OP2}$</td></tr><tr><td>Rate of change of labile organic phosphorus (OP<sub>1</sub>)</td><td>$\frac{dOP_1}{dt} = v_{Sub} \cdot I_{OP1} - OP_1Accumulation - OP_1Degraded - W_{OP1}$</td></tr><tr><td>Rate of change of refractory organic phosphorus (OP<sub>2</sub>)</td><td>$\frac{dOP_2}{dt} = v_{Sub} \cdot I_{OP2} - OP_2Accumulation - OP_2Degraded - W_{OP2}$</td></tr><tr><td colspan="2"><strong>Carbon cycle – biomass component</strong></td></tr><tr><td>Growth of subglacial autotrophs (A<sub>1</sub>)</td><td>$U_{A1} = A_1 \cdot T_f \cdot d \cdot p_{sub} \cdot I_{maxA} \cdot \left( \frac{DIN}{DIN+(K_n \cdot k_{Sub})} \right) \cdot \left( \frac{DIP}{DIP+(K_p \cdot k_{Sub})} \right)$</td></tr><tr><td>Growth of soil autotrophs (A<sub>2</sub>)</td><td>$U_{A2} = A_2 \cdot T_f \cdot d \cdot I_{maxA} \cdot \left( \frac{PAR}{PAR+K_L} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right)$</td></tr><tr><td>Growth of N-fixing soil autotrophs (A<sub>3</sub>)</td><td>$U_{A3} = U_{A3\_N2} + U_{A3\_DIN}$</td></tr><tr><td>Growth of N-fixing soil autotrophs (A<sub>3</sub>) with nitrogen fixation</td><td>$U_{A3\_N2} = A_3 \cdot T_f \cdot d \cdot n_f \cdot I_{maxA3} \cdot \left( \frac{PAR}{PAR+K_L} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (1 - m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil autotrophs (A<sub>3</sub>) with DIN</td><td>$U_{A3\_DIN} = A_3 \cdot T_f \cdot d \cdot I_{maxA3} \cdot \left( \frac{PAR}{PAR+K_L} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (m_{DIN})$</td></tr><tr><td>Growth of subglacial heterotrophs (H<sub>1</sub>) from labile substrate</td><td>$U_{H1L} = H_1 \cdot T_f \cdot d \cdot p_{sub} \cdot I_{maxH} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+(K_s \cdot k_{Sub})} \right) \cdot \left( \frac{DIN}{DIN+(K_n \cdot k_{Sub})} \right) \cdot \left( \frac{DIP}{DIP+(K_p \cdot k_{Sub})} \right)$</td></tr><tr><td>Growth of soil heterotrophs (H<sub>2</sub>) from labile substrate</td><td>$U_{H2L} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right)$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from labile substrate and nitrogen fixation</td><td>$U_{H3L\_N2} = H_3 \cdot T_f \cdot d \cdot n_f \cdot I_{maxH3} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+K_s} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (1 - m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from labile substrate and DIN</td><td>$U_{H3L\_DIN} = H_3 \cdot T_f \cdot d \cdot I_{maxH3} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (m_{DIN})$</td></tr><tr><td>Growth of subglacial heterotrophs (H<sub>1</sub>) from refractory substrate</td><td>$U_{H1R} = H_1 \cdot T_f \cdot d \cdot I_{maxH} \cdot p_{sub} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+(K_s \cdot k_{Sub})} \right) \cdot \left( \frac{DIN}{DIN+(K_n \cdot k_{Sub})} \right) \cdot \left( \frac{DIP}{DIP+(K_p \cdot k_{Sub})} \right)$</td></tr><tr><td>Growth of soil heterotrophs (H<sub>2</sub>) from refractory substrate</td><td>$U_{H2R} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right)$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from refractory substrate and nitrogen fixation</td><td>$U_{H3R\_N2} = H_3 \cdot T_f \cdot d \cdot n_f \cdot I_{maxH3} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+K_s} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (1 - m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from refractory substrate and DIN</td><td>$U_{H3R\_DIN} = H_3 \cdot T_f \cdot d \cdot I_{maxH3} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from labile substrate</td><td>$U_{H3L} = U_{H3L\_N2} + U_{H3L\_DIN}$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from refractory substrate</td><td>$U_{H3R} = U_{H3R\_N2} + U_{H3R\_DIN}$</td></tr></tbody></table>
<table><caption>**Table 4.** SHIMMER model formulation.</caption><thead><tr><th colspan="2">Fundamental balance equations</th></tr></thead><tbody><tr><td>Rate of change of autotrophic biomass (A<sub>1-3</sub>)</td><td>$\frac{dA_i}{dt} = U_{Ai} - G_{Ai} - X_{Ai}$</td></tr><tr><td>Rate of change of heterotrophic biomass (H<sub>1-3</sub>)</td><td>$\frac{dH_i}{dt} = U_{Hi} - G_{Hi} - X_{Hi}$</td></tr><tr><td>Rate of change of labile carbon substrate (S<sub>1</sub>)</td><td>$\frac{dS_1}{dt} = v_{Sub} \cdot I_{S1} + q \cdot G_{Ai} + q \cdot G_{Hi} + X_A + X_H - U_{S1Hi} - W_{S1}$</td></tr><tr><td>Rate of change of refractory carbon substrate (S<sub>2</sub>)</td><td>$\frac{dS_2}{dt} = v_{CS2} \cdot I_{S2} + (1-q) \cdot G_{Ai} + (1-q) \cdot G_{Hi} - U_{S2Hi} - W_{S2}$</td></tr><tr><td>Rate of change of dissolved inorganic nitrogen (DIN)</td><td>$\frac{dDIN}{dt} = v_{DIN} \cdot I_{DIN} - DIN_{Consumed} + DIN_{Released} - W_{DIN}$</td></tr><tr><td>Rate of change of dissolved inorganic phosphorus (DIP)</td><td>$\frac{dDIP}{dt} = v_{DIP} \cdot I_{DIP} - DIP_{Consumed} + DIP_{Released} - W_{DIP}$</td></tr><tr><td>Rate of change of labile organic nitrogen (ON<sub>1</sub>)</td><td>$\frac{dON_1}{dt} = v_{Sub} \cdot I_{ON1} - ON_1Accumulation - ON_1Degraded - W_{ON1}$</td></tr><tr><td>Rate of change of refractory organic nitrogen (ON<sub>2</sub>)</td><td>$\frac{dON_2}{dt} = v_{Sub} \cdot I_{OP2} - OP_2Accumulation - OP_2Degraded - W_{OP2}$</td></tr><tr><td>Rate of change of labile organic phosphorus (OP<sub>1</sub>)</td><td>$\frac{dOP_1}{dt} = v_{Sub} \cdot I_{OP1} - OP_1Accumulation - OP_1Degraded - W_{OP1}$</td></tr><tr><td>Rate of change of refractory organic phosphorus (OP<sub>2</sub>)</td><td>$\frac{dOP_2}{dt} = v_{Sub} \cdot I_{OP2} - OP_2Accumulation - OP_2Degraded - W_{OP2}$</td></tr><tr><td colspan="2">Carbon cycle – biomass component</td></tr><tr><td>Growth of subglacial autotrophs (A<sub>1</sub>)</td><td>$U_{A1} = A_1 \cdot T_f \cdot d \cdot p_{sub} \cdot I_{maxA} \cdot \left( \frac{DIN}{DIN+(K_n \cdot k_{Sub})} \right) \cdot \left( \frac{DIP}{DIP+(K_p \cdot k_{Sub})} \right)$</td></tr><tr><td>Growth of soil autotrophs (A<sub>2</sub>)</td><td>$U_{A2} = A_2 \cdot T_f \cdot d \cdot I_{maxA} \cdot \left( \frac{PAR}{PAR+K_L} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right)$</td></tr><tr><td>Growth of N-fixing soil autotrophs (A<sub>3</sub>)</td><td>$U_{A3} = U_{A3\_N2} + U_{A3\_DIN}$</td></tr><tr><td>Growth of N-fixing soil autotrophs (A<sub>3</sub>) with nitrogen fixation</td><td>$U_{A3\_N2} = A_3 \cdot T_f \cdot d \cdot n_f \cdot I_{maxA3} \cdot \left( \frac{PAR}{PAR+K_L} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (1 - m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil autotrophs (A<sub>3</sub>) with DIN</td><td>$U_{A3\_DIN} = A_3 \cdot T_f \cdot d \cdot I_{maxA3} \cdot \left( \frac{PAR}{PAR+K_L} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (m_{DIN})$</td></tr><tr><td>Growth of subglacial heterotrophs (H<sub>1</sub>) from labile substrate</td><td>$U_{H1L} = H_1 \cdot T_f \cdot d \cdot p_{sub} \cdot I_{maxH} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+(K_s \cdot k_{Sub})} \right) \cdot \left( \frac{DIN}{DIN+(K_n \cdot k_{Sub})} \right) \cdot \left( \frac{DIP}{DIP+(K_p \cdot k_{Sub})} \right)$</td></tr><tr><td>Growth of soil heterotrophs (H<sub>2</sub>) from labile substrate</td><td>$U_{H2L} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right)$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from labile substrate and nitrogen fixation</td><td>$U_{H3L\_N2} = H_3 \cdot T_f \cdot d \cdot n_f \cdot I_{maxH3} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+K_s} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (1 - m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from labile substrate and DIN</td><td>$U_{H3L\_DIN} = H_3 \cdot T_f \cdot d \cdot I_{maxH3} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (m_{DIN})$</td></tr><tr><td>Growth of subglacial heterotrophs (H<sub>1</sub>) from refractory substrate</td><td>$U_{H1R} = H_1 \cdot T_f \cdot d \cdot I_{maxH} \cdot p_{sub} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+(K_s \cdot k_{Sub})} \right) \cdot \left( \frac{DIN}{DIN+(K_n \cdot k_{Sub})} \right) \cdot \left( \frac{DIP}{DIP+(K_p \cdot k_{Sub})} \right)$</td></tr><tr><td>Growth of soil heterotrophs (H<sub>2</sub>) from refractory substrate</td><td>$U_{H2R} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right)$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from refractory substrate and nitrogen fixation</td><td>$U_{H3R\_N2} = H_3 \cdot T_f \cdot d \cdot n_f \cdot I_{maxH3} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+K_s} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (1 - m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from refractory substrate and DIN</td><td>$U_{H3R\_DIN} = H_3 \cdot T_f \cdot d \cdot I_{maxH3} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2+K_s} \right) \cdot \left( \frac{DIN}{DIN+K_n} \right) \cdot \left( \frac{DIP}{DIP+K_p} \right) \cdot (m_{DIN})$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from labile substrate</td><td>$U_{H3L} = U_{H3L\_N2} + U_{H3L\_DIN}$</td></tr><tr><td>Growth of N-fixing soil heterotrophs (H<sub>3</sub>) from refractory substrate</td><td>$U_{H3R} = U_{H3R\_N2} + U_{H3R\_DIN}$</td></tr></tbody></table>
<table><caption>**Table 4. Continued.**</caption><tbody><tr><td colspan="2"><strong>Fundamental balance equations</strong></td></tr><tr><td>Total growth of heterotrophs (H<sub>i</sub>)</td><td>U<sub>Hi</sub> = U<sub>HiL</sub> + U<sub>HiR</sub></td></tr><tr><td>Loss of autotrophic biomass</td><td>G<sub>Ai</sub> = T<sub>f</sub> &middot; d &middot; &alpha;<sub>A</sub> &middot; A<sub>i</sub><sup>2</sup></td></tr><tr><td>Loss of heterotrophic biomass</td><td>G<sub>Hi</sub> = T<sub>f</sub> &middot; d &middot; &alpha;<sub>H</sub> &middot; H<sub>i</sub><sup>2</sup></td></tr><tr><td>Total deaths that become labile</td><td>G<sub>L</sub> = q &middot; (G<sub>A1</sub> + G<sub>A2</sub> + G<sub>A3</sub> + G<sub>H1</sub> + G<sub>H2</sub> + G<sub>H3</sub>)</td></tr><tr><td>Total deaths that become refractory</td><td>G<sub>R</sub> = (1 &minus; q) &middot; (G<sub>A1</sub> + G<sub>A2</sub> + G<sub>A3</sub> + G<sub>H1</sub> + G<sub>H2</sub> + G<sub>H3</sub>)</td></tr><tr><td>Rate of exudate and EPS production by autotrophs</td><td>X<sub>Ai</sub> = ex<sub>A</sub> &middot; U<sub>Ai</sub></td></tr><tr><td>Rate of exudate and EPS production by heterotrophs</td><td>X<sub>Hi</sub> = ex<sub>H</sub> &middot; U<sub>Hi</sub></td></tr><tr><td>Total exudate production</td><td>X<sub>Total</sub> = (X<sub>A1</sub> + X<sub>A2</sub> + X<sub>A3</sub> + X<sub>H1</sub> + X<sub>H2</sub> + X<sub>H3</sub>)</td></tr><tr><td colspan="2"><strong>Carbon cycle – substrate component</strong></td></tr><tr><td>Consumption of labile substrate</td><td>S<sub>1Consumed</sub> = (1/Y<sub>H</sub>) &middot; U<sub>H1L</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H2L</sub> + (1/(n<sub>f</sub>&middot;Y<sub>H</sub>)) &middot; U<sub>H3L_N2</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H3L_DIN</sub></td></tr><tr><td>Consumption of refractory substrate</td><td>S<sub>2Consumed</sub> = (1/Y<sub>H</sub>) &middot; U<sub>H1R</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H2R</sub> + (1/(n<sub>f</sub>&middot;Y<sub>H</sub>)) &middot; U<sub>H3R_N2</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H3R_DIN</sub></td></tr><tr><td>Leaching</td><td>W<sub>X</sub> = g<sub>X</sub> &middot; X</td></tr><tr><td colspan="2"><strong>Nitrogen cycle</strong></td></tr><tr><td>Labile organic nitrogen degraded</td><td>ON<sub>1Degraded</sub> = (ON<sub>1</sub>/S<sub>1</sub>) &middot; (S<sub>1Consumed</sub>)</td></tr><tr><td>Refractory organic nitrogen degraded</td><td>ON<sub>2Degraded</sub> = (ON<sub>2</sub>/S<sub>2</sub>) &middot; (S<sub>2Consumed</sub>)</td></tr><tr><td>Labile organic nitrogen accumulated</td><td>ON<sub>1Accumulation</sub> = NC &middot; (G<sub>L</sub> + X<sub>T</sub>)</td></tr><tr><td>Refractory organic nitrogen accumulated</td><td>ON<sub>2Accumulation</sub> = NC &middot; G<sub>R</sub></td></tr><tr><td>DIN consumed</td><td>DIN<sub>Consumed</sub> = NC &middot; (U<sub>A1</sub> + U<sub>A2</sub> + U<sub>A3_DIN</sub> + U<sub>H1</sub> + U<sub>H2</sub> + U<sub>H3L_DIN</sub> + U<sub>H3R_DIN</sub>)</td></tr><tr><td>DIN released</td><td>DIN<sub>Released</sub> = ON<sub>1Degraded</sub> + ON<sub>2Degraded</sub></td></tr><tr><td colspan="2"><strong>Phosphorus cycle</strong></td></tr><tr><td>Labile organic phosphorus degraded</td><td>OP<sub>1Degraded</sub> = (OP<sub>1</sub>/S<sub>1</sub>) &middot; (S<sub>1Consumed</sub>)</td></tr><tr><td>Refractory organic phosphorus degraded</td><td>OP<sub>2Degraded</sub> = (OP<sub>2</sub>/S<sub>2</sub>) &middot; (S<sub>2Consumed</sub>)</td></tr><tr><td>Labile organic phosphorus accumulated</td><td>OP<sub>1Accumulation</sub> = PC &middot; (G<sub>L</sub> + X<sub>T</sub>)</td></tr><tr><td>Refractory organic phosphorus accumulated</td><td>OP<sub>2Accumulation</sub> = PC &middot; G<sub>R</sub></td></tr><tr><td>DIP consumed</td><td>DIP<sub>Consumed</sub> = PC &middot; (U<sub>A1</sub> + U<sub>A2</sub> + U<sub>A3</sub> + U<sub>H1</sub> + U<sub>H2</sub> + U<sub>H3</sub>)</td></tr><tr><td>DIP released</td><td>DIP<sub>Released</sub> = OP<sub>1Degraded</sub> + OP<sub>2Degraded</sub></td></tr><tr><td colspan="2"><strong>Environmental and scaling equations</strong></td></tr><tr><td>Temperature factor response</td><td>T<sub>f</sub> = exp((T &minus; T<sub>ref</sub>) / 10) &middot; log<sub>e</sub>(Q<sub>10</sub>)</td></tr><tr><td>Monod expression for nitrogen fixation inhibition in the presence of DIN</td><td>if DIN &le; DIN<sub>t</sub><br>m<sub>DIN</sub> = 0<br>else m<sub>DIN</sub> = ((DIN &minus; DIN<sub>t</sub>) / (DIN &minus; DIN<sub>t</sub>) + K<sub>N2</sub>)</td></tr></tbody></table>

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}$).

The nitrogen-fixing groups ($A_3$ and $H_3$) have the ability to fix atmospheric nitrogen gas under nitrogen-limiting conditions, or assimilate DIN. Bacteria in environmental samples

do not usually fix nitrogen in the presence of available DIN sources (Bottomley and Myrold, 2007; Holl and Montoya, 2005) because nitrogen fixation is an energetically expensive process (Liu et al., 2011; Cannell and Thornley, 2000; Phillips, 1980). However, even at very high concentrations of nitrate, there may not be total inhibition of nitrogen fixation, and simultaneous nitrogen fixation and nitrate assimilation may occur (Holl and Montoya, 2005). The inhibition of nitrogen fixation with DIN in SHIMMER thus follows Monod
<table><caption>**Table 4. Continued.**</caption><tbody><tr><td colspan="2"><strong>Fundamental balance equations</strong></td></tr><tr><td>Total growth of heterotrophs (H<sub>i</sub>)</td><td>U<sub>Hi</sub> = U<sub>HiL</sub> + U<sub>HiR</sub></td></tr><tr><td>Loss of autotrophic biomass</td><td>G<sub>Ai</sub> = T<sub>f</sub> &middot; d &middot; &alpha;<sub>A</sub> &middot; A<sub>i</sub><sup>2</sup></td></tr><tr><td>Loss of heterotrophic biomass</td><td>G<sub>Hi</sub> = T<sub>f</sub> &middot; d &middot; &alpha;<sub>H</sub> &middot; H<sub>i</sub><sup>2</sup></td></tr><tr><td>Total deaths that become labile</td><td>G<sub>L</sub> = q &middot; (G<sub>A1</sub> + G<sub>A2</sub> + G<sub>A3</sub> + G<sub>H1</sub> + G<sub>H2</sub> + G<sub>H3</sub>)</td></tr><tr><td>Total deaths that become refractory</td><td>G<sub>R</sub> = (1 &minus; q) &middot; (G<sub>A1</sub> + G<sub>A2</sub> + G<sub>A3</sub> + G<sub>H1</sub> + G<sub>H2</sub> + G<sub>H3</sub>)</td></tr><tr><td>Rate of exudate and EPS production by autotrophs</td><td>X<sub>Ai</sub> = ex<sub>A</sub> &middot; U<sub>Ai</sub></td></tr><tr><td>Rate of exudate and EPS production by heterotrophs</td><td>X<sub>Hi</sub> = ex<sub>H</sub> &middot; U<sub>Hi</sub></td></tr><tr><td>Total exudate production</td><td>X<sub>Total</sub> = (X<sub>A1</sub> + X<sub>A2</sub> + X<sub>A3</sub> + X<sub>H1</sub> + X<sub>H2</sub> + X<sub>H3</sub>)</td></tr><tr><td colspan="2"><strong>Carbon cycle – substrate component</strong></td></tr><tr><td>Consumption of labile substrate</td><td>S<sub>1Consumed</sub> = (1/Y<sub>H</sub>) &middot; U<sub>H1L</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H2L</sub> + (1/(n<sub>f</sub>&middot;Y<sub>H</sub>)) &middot; U<sub>H3L_N2</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H3L_DIN</sub></td></tr><tr><td>Consumption of refractory substrate</td><td>S<sub>2Consumed</sub> = (1/Y<sub>H</sub>) &middot; U<sub>H1R</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H2R</sub> + (1/(n<sub>f</sub>&middot;Y<sub>H</sub>)) &middot; U<sub>H3R_N2</sub> + (1/Y<sub>H</sub>) &middot; U<sub>H3R_DIN</sub></td></tr><tr><td>Leaching</td><td>W<sub>X</sub> = g<sub>X</sub> &middot; X</td></tr><tr><td colspan="2"><strong>Nitrogen cycle</strong></td></tr><tr><td>Labile organic nitrogen degraded</td><td>ON<sub>1Degraded</sub> = (ON<sub>1</sub>/S<sub>1</sub>) &middot; (S<sub>1Consumed</sub>)</td></tr><tr><td>Refractory organic nitrogen degraded</td><td>ON<sub>2Degraded</sub> = (ON<sub>2</sub>/S<sub>2</sub>) &middot; (S<sub>2Consumed</sub>)</td></tr><tr><td>Labile organic nitrogen accumulated</td><td>ON<sub>1Accumulation</sub> = NC &middot; (G<sub>L</sub> + X<sub>T</sub>)</td></tr><tr><td>Refractory organic nitrogen accumulated</td><td>ON<sub>2Accumulation</sub> = NC &middot; G<sub>R</sub></td></tr><tr><td>DIN consumed</td><td>DIN<sub>Consumed</sub> = NC &middot; (U<sub>A1</sub> + U<sub>A2</sub> + U<sub>A3_DIN</sub> + U<sub>H1</sub> + U<sub>H2</sub> + U<sub>H3L_DIN</sub> + U<sub>H3R_DIN</sub>)</td></tr><tr><td>DIN released</td><td>DIN<sub>Released</sub> = ON<sub>1Degraded</sub> + ON<sub>2Degraded</sub></td></tr><tr><td colspan="2"><strong>Phosphorus cycle</strong></td></tr><tr><td>Labile organic phosphorus degraded</td><td>OP<sub>1Degraded</sub> = (OP<sub>1</sub>/S<sub>1</sub>) &middot; (S<sub>1Consumed</sub>)</td></tr><tr><td>Refractory organic phosphorus degraded</td><td>OP<sub>2Degraded</sub> = (OP<sub>2</sub>/S<sub>2</sub>) &middot; (S<sub>2Consumed</sub>)</td></tr><tr><td>Labile organic phosphorus accumulated</td><td>OP<sub>1Accumulation</sub> = PC &middot; (G<sub>L</sub> + X<sub>T</sub>)</td></tr><tr><td>Refractory organic phosphorus accumulated</td><td>OP<sub>2Accumulation</sub> = PC &middot; G<sub>R</sub></td></tr><tr><td>DIP consumed</td><td>DIP<sub>Consumed</sub> = PC &middot; (U<sub>A1</sub> + U<sub>A2</sub> + U<sub>A3</sub> + U<sub>H1</sub> + U<sub>H2</sub> + U<sub>H3</sub>)</td></tr><tr><td>DIP released</td><td>DIP<sub>Released</sub> = OP<sub>1Degraded</sub> + OP<sub>2Degraded</sub></td></tr><tr><td colspan="2"><strong>Environmental and scaling equations</strong></td></tr><tr><td>Temperature factor response</td><td>T<sub>f</sub> = exp((T &minus; T<sub>ref</sub>) / 10) log<sub>e</sub>(Q<sub>10</sub>)</td></tr><tr><td>Monod expression for nitrogen fixation inhibition in the presence of DIN</td><td>if DIN &le; DIN<sub>t</sub><br>m<sub>DIN</sub> = 0<br>else m<sub>DIN</sub> = ((DIN &minus; DIN<sub>t</sub>) / (DIN &minus; DIN<sub>t</sub>) + K<sub>N2</sub>)</td></tr></tbody></table>

species (A<sub>1</sub> and H<sub>1</sub>) are represented in SHIMMER as more energy conserving and adapted to harsh environmental conditions. Accordingly, their maximum growth rate (I<sub>max</sub>) is reduced (by a factor p<sub>Sub</sub>), and their nutrient limitation response (expressed by the half-saturation constants K<sub>S</sub>, K<sub>N</sub> and K<sub>P</sub>) is reduced (by a factor k<sub>Sub</sub>).

The nitrogen-fixing groups (A<sub>3</sub> and H<sub>3</sub>) have the ability to fix atmospheric nitrogen gas under nitrogen-limiting conditions, or assimilate DIN. Bacteria in environmental samples

do not usually fix nitrogen in the presence of available DIN sources (Bottomley and Myrold, 2007; Holl and Montoya, 2005) because nitrogen fixation is an energetically expensive process (Liu et al., 2011; Cannell and Thornley, 2000; Phillips, 1980). However, even at very high concentrations of nitrate, there may not be total inhibition of nitrogen fixation, and simultaneous nitrogen fixation and nitrate assimilation may occur (Holl and Montoya, 2005). The inhibition of nitrogen fixation with DIN in SHIMMER thus follows Monod
<table><caption>Table 5. Parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Units</th><th>Nominal value (reference)</th><th>Lower range (reference)</th><th>Upper range (reference)</th><th>Damma Glacier, Switzerland</th><th>Athabasca Glacier, Canada</th></tr></thead><tbody><tr><td>T<sub>ref</sub></td><td>Reference temperature for rates</td><td>°C</td><td>25 (Frey et al., 2010)</td><td>Fixed</td><td>Fixed</td><td>–</td><td>–</td></tr><tr><td>NC</td><td>C : N ratio (mass)</td><td>Unitless</td><td>0.141 (Bernasconi et al., 2011)</td><td>Fixed</td><td>Fixed</td><td>–</td><td>–</td></tr><tr><td>PC</td><td>C : P ratio (mass)</td><td>Unitless</td><td>0.083 (Bernasconi et al., 2011)</td><td>Fixed</td><td>Fixed</td><td>–</td><td>–</td></tr><tr><td>Q<sub>10</sub></td><td>Temperature sensitivity</td><td>Unitless</td><td>2.0 (Soetaert and Herman, 2009)</td><td>1.5</td><td>3.1 (Yoshitake et al., 2010)</td><td>2.2</td><td>2.2</td></tr><tr><td>α<sub>A</sub></td><td>Death rate (autotrophs)</td><td>d<sup>-1</sup></td><td>0.0120 (German et al., 2012)</td><td>0.0060 (Grant et al., 1993)</td><td>0.4800 (Scott et al., 1995)</td><td>–</td><td>0.0340</td></tr><tr><td>α<sub>H</sub></td><td>Death rate (heterotrophs)</td><td>d<sup>-1</sup></td><td>0.0120 (German et al., 2012)</td><td>0.0060 (Grant et al., 1993)</td><td>0.4800 (Scott et al., 1995)</td><td>–</td><td>0.0340</td></tr><tr><td>I<sub>maxA</sub></td><td>Maximum growth rate (autotrophs)</td><td>d<sup>-1</sup></td><td>1.21 (Frey et al., 2010)</td><td>0.30 (Mur et al., 1999)</td><td>1.40 (Mur et al., 1999)</td><td>1.40</td><td>1.4</td></tr><tr><td>I<sub>maxH</sub></td><td>Maximum growth rate (heterotrophs)</td><td>d<sup>-1</sup></td><td>1.21 (Frey et al., 2010)</td><td>0.24 (Ingwersen et al., 2008)</td><td>4.80 (Darrah, 1991; Blagodatsky and Richter, 1998)</td><td>1.24</td><td>1.24</td></tr><tr><td>ex<sub>A</sub></td><td>Exudates and EPS production (autotrophs)</td><td>Unitless</td><td>0.014 (Allison, 2005)</td><td>0.007 (Allison, 2005)</td><td>0.021 (Allison, 2005)</td><td>–</td><td>–</td></tr><tr><td>ex<sub>H</sub></td><td>Exudates and EPS production (heterotrophs)</td><td>Unitless</td><td>0.014 (Allison, 2005)</td><td>0.007 (Allison, 2005)</td><td>0.021 (Allison, 2005)</td><td>–</td><td>–</td></tr><tr><td>p<sub>Sub</sub></td><td>Slow down of subglacial microbial growth rate</td><td>Unitless</td><td>0.2</td><td>0.1</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>k<sub>Sub</sub></td><td>Lower half-saturation constants (K<sub>S</sub>, K<sub>N</sub> and K<sub>P</sub>) for subglacial microbes</td><td>Unitless</td><td>0.8</td><td>0.1</td><td>1.0</td><td>–</td><td>–</td></tr><tr><td>K<sub>L</sub></td><td>Light half-saturation constant for autotrophs (A<sub>2</sub> and A<sub>3</sub>)</td><td>W m<sup>-2</sup> (PAR)</td><td>1.85 (Van Liere and Walsby, 1982)</td><td>0.70 (Van Liere and Walsby, 1982)</td><td>3.00 (Van Liere and Walsby, 1982)</td><td>–</td><td>–</td></tr><tr><td>K<sub>S</sub></td><td>Substrate half-saturation constant for heterotrophs</td><td>µg g<sup>-1</sup></td><td>349 (Vandewerf and Verstraete, 1987b)</td><td>50 (Darrah, 1991)</td><td>1000 (Knapp et al., 1983)</td><td>–</td><td>–</td></tr><tr><td>K<sub>N</sub></td><td>DIN half-saturation constant</td><td>µg g<sup>-1</sup></td><td>49.209 (stoichiometric)</td><td>7.050 (stoichiometric)</td><td>141.000 (stoichiometric)</td><td>–</td><td>–</td></tr><tr><td>K<sub>P</sub></td><td>DIP half-saturation constant</td><td>µg g<sup>-1</sup></td><td>28.967 (stoichiometric)</td><td>4.150 (stoichiometric)</td><td>83.000 (stoichiometric)</td><td>–</td><td>–</td></tr><tr><td>n<sub>f</sub></td><td>Downscaling of Y and I<sub>max</sub> when fixing nitrogen</td><td>Unitless</td><td>0.50 (Bottomley and Myrold, 2007)</td><td>0.10 (LaRoche and Breitbarth, 2005; Breitbarth et al., 2008; Goebel et al., 2008)</td><td>0.80</td><td>–</td><td>–</td></tr><tr><td>K<sub>N2</sub></td><td>Nitrogen fixation inhibition</td><td>µg g<sup>-1</sup></td><td>393.672 (Holl and Montoya, 2005; Rabouille et al., 2006)</td><td>56.4 (Holl and Montoya, 2005; Rabouille et al., 2006)</td><td>1128 (Holl and Montoya, 2005; Rabouille et al., 2006)</td><td>–</td><td>–</td></tr><tr><td>DIN<sub>f</sub></td><td>Threshold value of DIN for nitrogen fixation inhibition</td><td>µg g<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>–</td><td>–</td></tr><tr><td>q</td><td>Proportion of necromass that becomes labile (S<sub>1</sub>)</td><td>Unitless</td><td>0.3</td><td>0.1</td><td>0.5</td><td>–</td><td>–</td></tr><tr><td>J<sub>S1</sub></td><td>Bioavailability (preference) of S<sub>1</sub></td><td>Unitless</td><td>0.50</td><td>0.50</td><td>1.00</td><td>0.68</td><td>0.50</td></tr><tr><td>J<sub>S2</sub></td><td>Bioavailability (preference) of S<sub>2</sub></td><td>Unitless</td><td>0.10</td><td>0.10</td><td>0.50</td><td>0.15</td><td>0.10</td></tr><tr><td>g<sub>Sub</sub></td><td>Leaching of substrate</td><td>d<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>–</td><td>–</td></tr><tr><td>g<sub>DIN</sub></td><td>Leaching of DIN</td><td>d<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>–</td><td>–</td></tr><tr><td>g<sub>DIP</sub></td><td>Leaching of DIP</td><td>d<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>–</td><td>–</td></tr><tr><td>Y<sub>A</sub></td><td>Growth efficiency of autotrophs</td><td>g C (g C consumed)<sup>-1</sup></td><td>0.200 (Scott et al., 1995)</td><td>0.100 (Foereid and Yearsley, 2004)</td><td>0.848 (Blagodatsky and Richter, 1998)</td><td>–</td><td>–</td></tr><tr><td>Y<sub>H</sub></td><td>Growth efficiency of heterotrophs</td><td>g C (g C consumed)<sup>-1</sup></td><td>0.200 (Scott et al., 1995)</td><td>0.134 (German et al., 2012)</td><td>0.848 (Blagodatsky and Richter, 1998)</td><td>–</td><td>–</td></tr></tbody></table>
<table><caption>Table 5. Parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Units</th><th>Nominal value (reference)</th><th>Lower range (reference)</th><th>Upper range (reference)</th><th>Damma Glacier, Switzerland</th><th>Athabasca Glacier, Canada</th></tr></thead><tbody><tr><td>T<sub>ref</sub></td><td>Reference temperature for rates</td><td>°C</td><td>25 (Frey et al., 2010)</td><td>Fixed</td><td>Fixed</td><td>—</td><td>—</td></tr><tr><td>NC</td><td>C : N ratio (mass)</td><td>Unitless</td><td>0.141 (Bernasconi et al., 2011)</td><td>Fixed</td><td>Fixed</td><td>—</td><td>—</td></tr><tr><td>PC</td><td>C : P ratio (mass)</td><td>Unitless</td><td>0.083 (Bernasconi et al., 2011)</td><td>Fixed</td><td>Fixed</td><td>—</td><td>—</td></tr><tr><td>Q<sub>10</sub></td><td>Temperature sensitivity</td><td>Unitless</td><td>2.0 (Soetaert and Herman, 2009)</td><td>1.5</td><td>3.1 (Yoshitake et al., 2010)</td><td>2.2</td><td>2.2</td></tr><tr><td>α<sub>A</sub></td><td>Death rate (autotrophs)</td><td>d<sup>-1</sup></td><td>0.0120 (German et al., 2012)</td><td>0.0060 (Grant et al., 1993)</td><td>0.4800 (Scott et al., 1995)</td><td>—</td><td>0.0340</td></tr><tr><td>α<sub>H</sub></td><td>Death rate (heterotrophs)</td><td>d<sup>-1</sup></td><td>0.0120 (German et al., 2012)</td><td>0.0060 (Grant et al., 1993)</td><td>0.4800 (Scott et al., 1995)</td><td>—</td><td>0.0340</td></tr><tr><td>I<sub>maxA</sub></td><td>Maximum growth rate (autotrophs)</td><td>d<sup>-1</sup></td><td>1.21 (Frey et al., 2010)</td><td>0.30 (Mur et al., 1999)</td><td>1.40 (Mur et al., 1999)</td><td>1.40</td><td>1.4</td></tr><tr><td>I<sub>maxH</sub></td><td>Maximum growth rate (heterotrophs)</td><td>d<sup>-1</sup></td><td>1.21 (Frey et al., 2010)</td><td>0.24 (Ingwersen et al., 2008)</td><td>4.80 (Darrah, 1991; Blagodatsky and Richter, 1998)</td><td>1.24</td><td>1.24</td></tr><tr><td>ex<sub>A</sub></td><td>Exudates and EPS production (autotrophs)</td><td>Unitless</td><td>0.014 (Allison, 2005)</td><td>0.007 (Allison, 2005)</td><td>0.021 (Allison, 2005)</td><td>—</td><td>—</td></tr><tr><td>ex<sub>H</sub></td><td>Exudates and EPS production (heterotrophs)</td><td>Unitless</td><td>0.014 (Allison, 2005)</td><td>0.007 (Allison, 2005)</td><td>0.021 (Allison, 2005)</td><td>—</td><td>—</td></tr><tr><td>p<sub>Sub</sub></td><td>Slow down of subglacial microbial growth rate</td><td>Unitless</td><td>0.2</td><td>0.1</td><td>1.0</td><td>—</td><td>—</td></tr><tr><td>k<sub>Sub</sub></td><td>Lower half-saturation constants (K<sub>S</sub>, K<sub>N</sub> and K<sub>P</sub>) for subglacial microbes</td><td>Unitless</td><td>0.8</td><td>0.1</td><td>1.0</td><td>—</td><td>—</td></tr><tr><td>K<sub>L</sub></td><td>Light half-saturation constant for autotrophs (A<sub>2</sub> and A<sub>3</sub>)</td><td>W m<sup>-2</sup> (PAR)</td><td>1.85 (Van Liere and Walsby, 1982)</td><td>0.70 (Van Liere and Walsby, 1982)</td><td>3.00 (Van Liere and Walsby, 1982)</td><td>—</td><td>—</td></tr><tr><td>K<sub>S</sub></td><td>Substrate half-saturation constant for heterotrophs</td><td>µg g<sup>-1</sup></td><td>349 (Vandewerf and Verstraete, 1987b)</td><td>50 (Darrah, 1991)</td><td>1000 (Knapp et al., 1983)</td><td>—</td><td>—</td></tr><tr><td>K<sub>N</sub></td><td>DIN half-saturation constant</td><td>µg g<sup>-1</sup></td><td>49.209 (stoichiometric)</td><td>7.050 (stoichiometric)</td><td>141.000 (stoichiometric)</td><td>—</td><td>—</td></tr><tr><td>K<sub>P</sub></td><td>DIP half-saturation constant</td><td>µg g<sup>-1</sup></td><td>28.967 (stoichiometric)</td><td>4.150 (stoichiometric)</td><td>83.000 (stoichiometric)</td><td>—</td><td>—</td></tr><tr><td>n<sub>f</sub></td><td>Downscaling of Y and I<sub>max</sub> when fixing nitrogen</td><td>Unitless</td><td>0.50 (Bottomley and Myrold, 2007)</td><td>0.10 (LaRoche and Breitbarth, 2005; Breitbarth et al., 2008; Goebel et al., 2008)</td><td>0.80</td><td>—</td><td>—</td></tr><tr><td>K<sub>N2</sub></td><td>Nitrogen fixation inhibition</td><td>µg g<sup>-1</sup></td><td>393.672 (Holl and Montoya, 2005; Rabouille et al., 2006)</td><td>56.4 (Holl and Montoya, 2005; Rabouille et al., 2006)</td><td>1128 (Holl and Montoya, 2005; Rabouille et al., 2006)</td><td>—</td><td>—</td></tr><tr><td>DIN<sub>f</sub></td><td>Threshold value of DIN for nitrogen fixation inhibition</td><td>µg g<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>—</td><td>—</td></tr><tr><td>q</td><td>Proportion of necromass that becomes labile (S<sub>1</sub>)</td><td>Unitless</td><td>0.3</td><td>0.1</td><td>0.5</td><td>—</td><td>—</td></tr><tr><td>J<sub>S1</sub></td><td>Bioavailability (preference) of S<sub>1</sub></td><td>Unitless</td><td>0.50</td><td>0.50</td><td>1.00</td><td>0.68</td><td>0.50</td></tr><tr><td>J<sub>S2</sub></td><td>Bioavailability (preference) of S<sub>2</sub></td><td>Unitless</td><td>0.10</td><td>0.10</td><td>0.50</td><td>0.15</td><td>0.10</td></tr><tr><td>g<sub>Sub</sub></td><td>Leaching of substrate</td><td>d<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>—</td><td>—</td></tr><tr><td>g<sub>DIN</sub></td><td>Leaching of DIN</td><td>d<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>—</td><td>—</td></tr><tr><td>g<sub>DIP</sub></td><td>Leaching of DIP</td><td>d<sup>-1</sup></td><td>0</td><td>0</td><td>0</td><td>—</td><td>—</td></tr><tr><td>Y<sub>A</sub></td><td>Growth efficiency of autotrophs</td><td>g C (g C consumed)<sup>-1</sup></td><td>0.200 (Scott et al., 1995)</td><td>0.100 (Foereid and Yearsley, 2004)</td><td>0.848 (Blagodatsky and Richter, 1998)</td><td>—</td><td>—</td></tr><tr><td>Y<sub>H</sub></td><td>Growth efficiency of heterotrophs</td><td>g C (g C consumed)<sup>-1</sup></td><td>0.200 (Scott et al., 1995)</td><td>0.134 (German et al., 2012)</td><td>0.848 (Blagodatsky and Richter, 1998)</td><td>—</td><td>—</td></tr></tbody></table>
<table><caption>**Table 5. Continued.**</caption><thead><tr><th>Parameter</th><th>Description</th><th>Units</th><th>Nominal value (reference)</th><th>Lower range (reference)</th><th>Upper range (reference)</th><th>Damma Glacier, Switzerland</th><th>Athabasca Glacier, Canada</th></tr></thead><tbody><tr><td>d</td><td>Active fraction of microbial biomass</td><td>Unitless</td><td>0.285 (Wang et al., 2014)</td><td>0.100 (Blagodatsky and Richter, 1998)</td><td>0.580 (Blagodatsky and Richter, 1998)</td><td>0.285</td><td>0.100</td></tr><tr><td>v<sub>Sub</sub></td><td>Proportion of allochthonous substrate deposition retained</td><td>Unitless</td><td>0.60</td><td>0.00</td><td>1.00</td><td>—</td><td>0.05</td></tr><tr><td>v<sub>DIN</sub></td><td>Proportion of allochthonous DIN deposition retained</td><td>Unitless</td><td>0.01</td><td>0.00</td><td>1.00</td><td>—</td><td>—</td></tr><tr><td>v<sub>DIP</sub></td><td>Proportion of allochthonous DIP deposition retained</td><td>Unitless</td><td>0.01</td><td>0.00</td><td>1.00</td><td>—</td><td>—</td></tr></tbody></table>

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_{H3\_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 growth-limiting half-saturation expression for soil DIN ($K_n$) is discounted from $U_{A3\_N2}$ and $U_{H3\_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_{H2L}$) and refractory ($U_{H2R}$) organic matter is formulated in a similar way to Eqs. (3) and (4) but depends on the bioavailability of organic matter rather than light:

$$U_{H2L} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1 + K_s} \right) \cdot \left( \frac{\text{DIN}}{\text{DIN} + K_n} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + K_p} \right), \quad (5)$$

$$U_{H2R} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2 + K_s} \right) \cdot \left( \frac{\text{DIN}}{\text{DIN} + K_n} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + K_p} \right). \quad (6)$$

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_{Hi}$) 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_{Hi}$) are density dependent and are also sensitive to variations in soil temperature ($T$), the active fraction $d$ and adjustable rate parameters $\alpha_A$ (autotrophs) and $\alpha_H$ (heterotrophs):

$$G_{Ai} = T_f \cdot d \cdot \alpha_A \cdot A_i^2, \quad (7)$$

$$G_{Hi} = T_f \cdot d \cdot \alpha_H \cdot H_i^2. \quad (8)$$

It is assumed that losses ($G_{Ai}$ and $G_{Hi}$) form insoluble microbial necromass (organic matter), comprising of organic
<table><caption>**Table 5. Continued.**</caption><thead><tr><th>Parameter</th><th>Description</th><th>Units</th><th>Nominal value (reference)</th><th>Lower range (reference)</th><th>Upper range (reference)</th><th>Damma Glacier, Switzerland</th><th>Athabasca Glacier, Canada</th></tr></thead><tbody><tr><td>d</td><td>Active fraction of microbial biomass</td><td>Unitless</td><td>0.285 (Wang et al., 2014)</td><td>0.100 (Blagodatsky and Richter, 1998)</td><td>0.580 (Blagodatsky and Richter, 1998)</td><td>0.285</td><td>0.100</td></tr><tr><td>v<sub>Sub</sub></td><td>Proportion of allochthonous substrate deposition retained</td><td>Unitless</td><td>0.60</td><td>0.00</td><td>1.00</td><td>—</td><td>0.05</td></tr><tr><td>v<sub>DIN</sub></td><td>Proportion of allochthonous DIN deposition retained</td><td>Unitless</td><td>0.01</td><td>0.00</td><td>1.00</td><td>—</td><td>—</td></tr><tr><td>v<sub>DIP</sub></td><td>Proportion of allochthonous DIP deposition retained</td><td>Unitless</td><td>0.01</td><td>0.00</td><td>1.00</td><td>—</td><td>—</td></tr></tbody></table>

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_{H3\_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 growth-limiting half-saturation expression for soil DIN ($K_n$) is discounted from $U_{A3\_N2}$ and $U_{H3\_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_{H2L}$) and refractory ($U_{H2R}$) organic matter is formulated in a similar way to Eqs. (3) and (4) but depends on the bioavailability of organic matter rather than light:

$$U_{H2L} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S1} \cdot \frac{S_1}{S_1 + K_s} \right) \cdot \left( \frac{\text{DIN}}{\text{DIN} + K_n} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + K_p} \right), \quad (5)$$

$$U_{H2R} = H_2 \cdot T_f \cdot d \cdot I_{maxH} \cdot \left( J_{S2} \cdot \frac{S_2}{S_2 + K_s} \right) \cdot \left( \frac{\text{DIN}}{\text{DIN} + K_n} \right) \cdot \left( \frac{\text{DIP}}{\text{DIP} + K_p} \right). \quad (6)$$

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_{Hi}$) 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_{Hi}$) are density dependent and are also sensitive to variations in soil temperature ($T$), the active fraction $d$ and adjustable rate parameters $\alpha_A$ (autotrophs) and $\alpha_H$ (heterotrophs):

$$G_{Ai} = T_f \cdot d \cdot \alpha_A \cdot A_i^2, \quad (7)$$

$$G_{Hi} = T_f \cdot d \cdot \alpha_H \cdot H_i^2. \quad (8)$$

It is assumed that losses ($G_{Ai}$ and $G_{Hi}$) form insoluble microbial necromass (organic matter), comprising of organic
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_{Hi}$) (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):

$$X_{Ai} = ex_A \cdot U_{Ai}, \qquad (9)$$

$$X_{Hi} = ex_H \cdot U_{Hi}. \qquad (10)$$

Exudate is highly reactive; therefore these stocks directly contribute to the labile substrate pools ($S_1$, ON$_1$, OP$_1$).

### 2.2.2 Organic matter dynamics

Organic matter dynamics are described by the following equations:

$$\frac{dS_1}{dt} = v_{\text{Sub}} \cdot I_{S1} + q \cdot G_{Ai} + q \cdot G_{Hi} + X_A + X_H - U_{S1Hi} - W_{S1}, \qquad (11)$$

$$\frac{dS_2}{dt} = v_{\text{Sub}} \cdot I_{S2} + (1-q) \cdot G_{Ai} + (1-q) \cdot G_{Hi} - U_{S2Hi} - W_{S2}. \qquad (12)$$

Organic matter accumulates in the soil due to microbial loss ($G_{Ai}$ and $G_{Hi}$), exudate and EPS production ($X_{Ai}$ and $X_{Hi}$), and allochthonous external carbon substrate inputs ($I_{S1}$ and $I_{S2}$), and is depleted due to heterotrophic growth ($U_{S1Hi}$ and $U_{S2Hi}$) 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_{S1Hi}$ and $U_{S2Hi}$) 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$).

### 2.2.3 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

$$\frac{d\text{DIN}}{dt} = v_{\text{DIN}} \cdot I_{\text{DIN}} - \text{DIN}_{\text{Consumed}} + \text{DIN}_{\text{Released}} - W_{\text{DIN}}.$$

(13)

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).

## 3 Implementation

### 3.1 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 “Isoda” 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.

### 3.2 Forcing and initial conditions

The following external forcings drive and regulate the system's dynamics:

- PAR (wavelength of approximately 400–700 nm) ($W \text{ m}^{-2}$)

- snow depth (m)

- soil temperature (°C)

- allochthonous nutrient inputs ($\mu g \text{ g}^{-1} \text{ d}^{-1}$)

- nutrient leaching rate.

The presence of snow attenuates sunlight and inhibits PAR from reaching the soil surface. This is accounted for in pre-processing of forcing data according to the equation:

$$n = n_0 e^{-mx}, \qquad (14)$$

whereby $n$ is the irradiance ($W \text{ m}^{-2}$), $x$ is the snow depth (m) and $m$ is the extinction coefficient for snow ($\text{m}^{-1}$). The extinction coefficients for various types of snow can be measured and an estimate of $6 \text{ 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 Table 2, respectively. They are directly informed from field studies (Sect. 4.3.1 and 4.3.2).
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_{Hi}$) (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):

$$X_{Ai} = ex_A \cdot U_{Ai}, \qquad (9)$$

$$X_{Hi} = ex_H \cdot U_{Hi}. \qquad (10)$$

Exudate is highly reactive; therefore these stocks directly contribute to the labile substrate pools ($S_1$, ON$_1$, OP$_1$).

### 2.2.2 Organic matter dynamics

Organic matter dynamics are described by the following equations:

$$\frac{dS_1}{dt} = v_{\text{Sub}} \cdot I_{S1} + q \cdot G_{Ai} + q \cdot G_{Hi} + X_A + X_H - U_{S1Hi} - W_{S1}, \qquad (11)$$

$$\frac{dS_2}{dt} = v_{\text{Sub}} \cdot I_{S2} + (1-q) \cdot G_{Ai} + (1-q) \cdot G_{Hi} - U_{S2Hi} - W_{S2}. \qquad (12)$$

Organic matter accumulates in the soil due to microbial loss ($G_{Ai}$ and $G_{Hi}$), exudate and EPS production ($X_{Ai}$ and $X_{Hi}$), and allochthonous external carbon substrate inputs ($I_{S1}$ and $I_{S2}$), and is depleted due to heterotrophic growth ($U_{S1Hi}$ and $U_{S2Hi}$) 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_{S1Hi}$ and $U_{S2Hi}$) 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$).

### 2.2.3 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

$$\frac{d\text{DIN}}{dt} = v_{\text{DIN}} \cdot I_{\text{DIN}} - \text{DIN}_{\text{Consumed}} + \text{DIN}_{\text{Released}} - W_{\text{DIN}}.$$

(13)

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).

## 3 Implementation

### 3.1 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 “Isoda” 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.

### 3.2 Forcing and initial conditions

The following external forcings drive and regulate the system's dynamics:

- PAR (wavelength of approximately 400–700 nm) ($W \text{ m}^{-2}$)

- snow depth (m)

- soil temperature (°C)

- allochthonous nutrient inputs ($\mu g \text{ g}^{-1} \text{ d}^{-1}$)

- nutrient leaching rate.

The presence of snow attenuates sunlight and inhibits PAR from reaching the soil surface. This is accounted for in pre-processing of forcing data according to the equation:

$$n = n_0 e^{-mx}, \qquad (14)$$

whereby $n$ is the irradiance ($W \text{ m}^{-2}$), $x$ is the snow depth (m) and $m$ is the extinction coefficient for snow ($\text{m}^{-1}$). The extinction coefficients for various types of snow can be measured and an estimate of $6 \text{ 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 Table 2, respectively. They are directly informed from field studies (Sect. 4.3.1 and 4.3.2).
<figure><img src="image_6.png" /><figcaption>Figure 4. Annual meteorological forcings for SHIMMER for the Damma Glacier, Switzerland (red), and the Athabasca Glacier, Canada (black). Forcing data are provided by the WSL (Switzerland), Achuff and Coen (1980) and Staveley AAFC (Canada).</figcaption></figure>

# 4 Model runs

## 4.1 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).

## 4.2 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)

- 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).

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

$$ \lambda(X, p) = \frac{p}{X} \cdot \frac{\delta X}{\delta p}, \qquad (15) $$

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
<figure><img src="image_7.png" /><figcaption>Figure 4. Annual meteorological forcings for SHIMMER for the Damma Glacier, Switzerland (red), and the Athabasca Glacier, Canada (black). Forcing data are provided by the WSL (Switzerland), Achuff and Coen (1980) and Staveley AAFC (Canada).</figcaption></figure>

# 4 Model runs

## 4.1 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).

## 4.2 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)

- 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).

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

$$ \lambda(X, p) = \frac{p}{X} \cdot \frac{\delta X}{\delta p}, \qquad (15) $$

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
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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-

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Simulated response of autotrophic (A1-3) biomass (blue line) and heterotrophic (H1-3) biomass (red line) to variation in (a) Q10 and (b) ImaxH 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.</figcaption></figure>

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

$$ \emptyset = \frac{X_{\text{max}} - X_{\text{min}}}{X} \times 100, \qquad (16) $$
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

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-

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Simulated response of autotrophic (A1-3) biomass (blue line) and heterotrophic (H1-3) biomass (red line) to variation in (a) Q10 and (b) ImaxH 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.</figcaption></figure>

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

$$ \emptyset = \frac{X_{\max} - X_{\min}}{X} \times 100, \qquad (16) $$
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.

## 4.3 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.

### 4.3.1 Case study 1: Damma Glacier, Switzerland

Published data sets of the biogeochemical development of the Damma Glacier forefield in Canton Uri, Switzerland (46.6° N, 8.5° E), are used to test and validate the model, and explore detailed model behaviour (Bernasconi et al., 2011). Over the last 2 decades, plant and microbial succession at this site has been extensively studied. Comprehensive data sets have been collected as part of the BigLink project (Bernasconi et al., 2011), with further detailed studies on nutrient cycling (Brankatschk et al., 2011; Bernasconi et al., 2011; Guelland et al., 2013a; Göransson et al., 2014; Smittenberg et al., 2012; Tamburini et al., 2012), microbial community composition (Duc et al., 2009; Frey et al., 2013; Lazzaro et al., 2012; Meola et al., 2014; Sigler and Zeyer, 2002; Zumsteg et al., 2012), soil depths (Rime et al., 2014) and soil activity (Goransson et al., 2011; Guelland et al., 2013b; Zumsteg et al., 2011). Therefore, the site is highly appropriate

to gain insight into model behaviour and biogeochemical processes.

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 \mu g \text{ C g}^{-1}$ in recently exposed soils to around $30 000 \mu g \text{ 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):

- total microbial biomass ($A_1+A_2+A_3+H_1+H_2+H_3$): presented as $C_{mic}$;

- carbon substrate ($S_1+S_2$): calculated as $\text{TOC}-C_{mic}$;

- ON ($\text{ON}_1+\text{ON}_2$): calculated as $\text{TN}-N_{mic}$;

- available DIP: presented as $P_{\text{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 5-year 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. Least-square error calculation and minimisation of errors are done only on those data points. The remaining data points from the later stages of succession are used as a test to see if microbial abundance in older soils falls within a plausible range, given the scatter in the observational data. A secondary data
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.

## 4.3 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.

### 4.3.1 Case study 1: Damma Glacier, Switzerland

Published data sets of the biogeochemical development of the Damma Glacier forefield in Canton Uri, Switzerland (46.6° N, 8.5° E), are used to test and validate the model, and explore detailed model behaviour (Bernasconi et al., 2011). Over the last 2 decades, plant and microbial succession at this site has been extensively studied. Comprehensive data sets have been collected as part of the BigLink project (Bernasconi et al., 2011), with further detailed studies on nutrient cycling (Brankatschk et al., 2011; Bernasconi et al., 2011; Guelland et al., 2013a; Göransson et al., 2014; Smittenberg et al., 2012; Tamburini et al., 2012), microbial community composition (Duc et al., 2009; Frey et al., 2013; Lazzaro et al., 2012; Meola et al., 2014; Sigler and Zeyer, 2002; Zumsteg et al., 2012), soil depths (Rime et al., 2014) and soil activity (Goransson et al., 2011; Guelland et al., 2013b; Zumsteg et al., 2011). Therefore, the site is highly appropriate

to gain insight into model behaviour and biogeochemical processes.

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 \mu g \text{ C g}^{-1}$ in recently exposed soils to around $30 000 \mu g \text{ 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):

- total microbial biomass ($A_1+A_2+A_3+H_1+H_2+H_3$): presented as $C_{mic}$;

- carbon substrate ($S_1+S_2$): calculated as $\text{TOC}-C_{mic}$;

- ON ($\text{ON}_1+\text{ON}_2$): calculated as $\text{TN}-N_{mic}$;

- available DIP: presented as $P_{\text{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 5-year 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. Least-square error calculation and minimisation of errors are done only on those data points. The remaining data points from the later stages of succession are used as a test to see if microbial abundance in older soils falls within a plausible range, given the scatter in the observational data. A secondary data
set of DIP accumulation complements the BigLink data set (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⁻² s⁻¹) which are converted to PAR (W m⁻²) by a conversion factor (0.219). PAR (W m⁻²), 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).

Allochthonous inputs to the Damma Glacier forefield are estimated in Brankatschk et al. (2011) based on chemical analyses of the snowpack and model simulations:

- carbon substrate: 75 µg C cm⁻² yr⁻¹

- DIN: 80 µg N cm⁻² yr⁻¹

- ON: 6.3 µg N cm⁻² yr⁻¹.

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, 1g dry soil occupies a surface area of 0.869 cm² (Guelland et al., 2013b). Since substrate material and DIN is liberated when the snow-pack 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⁻² yr⁻¹ (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 parameters $v_{Sub}$ (for all substrate pools), $v_{DIN}$ (for DIN) and $v_{DIP}$ (for DIP).

<table><caption>**Table 6.** Allochthonous deposition.</caption><thead><tr><th>Nutrient species</th><th>Input (µg g<sup>-1</sup> d<sup>-1</sup>)</th></tr></thead><tbody><tr><td>$S_1$</td><td>1.955</td></tr><tr><td>$S_2$</td><td>4.562</td></tr><tr><td>DIN</td><td>6.952</td></tr><tr><td>DIP</td><td>1.738</td></tr><tr><td>ON<sub>1</sub></td><td>0.164</td></tr><tr><td>ON<sub>2</sub></td><td>0.383</td></tr><tr><td>OP<sub>1</sub></td><td>0.162</td></tr><tr><td>OP<sub>2</sub></td><td>0.378</td></tr></tbody></table>

### 4.3.2 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):

- Total microbial biomass ($A_1+A_2+A_3+H_1+H_2+H_3$);

- Carbon substrate ($S_1+S_2$): calculated as $C_{org}$ – microbial biomass.

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) and linearly interpolated to
set of DIP accumulation complements the BigLink data set (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⁻² s⁻¹) which are converted to PAR (W m⁻²) by a conversion factor (0.219). PAR (W m⁻²), 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).

Allochthonous inputs to the Damma Glacier forefield are estimated in Brankatschk et al. (2011) based on chemical analyses of the snowpack and model simulations:

- carbon substrate: 75 µg C cm⁻² yr⁻¹

- DIN: 80 µg N cm⁻² yr⁻¹

- ON: 6.3 µg N cm⁻² yr⁻¹.

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, 1g dry soil occupies a surface area of 0.869 cm² (Guelland et al., 2013b). Since substrate material and DIN is liberated when the snow-pack 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⁻² yr⁻¹ (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 parameters $v_{Sub}$ (for all substrate pools), $v_{DIN}$ (for DIN) and $v_{DIP}$ (for DIP).

<table><caption>**Table 6.** Allochthonous deposition.</caption><thead><tr><th>Nutrient species</th><th>Input (µg g<sup>-1</sup> d<sup>-1</sup>)</th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>1.955</td></tr><tr><td>S<sub>2</sub></td><td>4.562</td></tr><tr><td>DIN</td><td>6.952</td></tr><tr><td>DIP</td><td>1.738</td></tr><tr><td>ON<sub>1</sub></td><td>0.164</td></tr><tr><td>ON<sub>2</sub></td><td>0.383</td></tr><tr><td>OP<sub>1</sub></td><td>0.162</td></tr><tr><td>OP<sub>2</sub></td><td>0.378</td></tr></tbody></table>

### 4.3.2 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):

- Total microbial biomass ($A_1+A_2+A_3+H_1+H_2+H_3$);

- Carbon substrate ($S_1 + S_2$): calculated as $C_{org} - \text{microbial biomass}$.

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) and linearly interpolated to
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Model output of the Damma Glacier forefield for a 75 year time series with nominal parameter values (Table 4).</figcaption></figure>

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_{\text{Sub}}$, $v_{\text{DIN}}$ and $v_{\text{DIP}}$ are adjusted. As with case study 1, optimisation is carried out based on the results of the sensitivity study, and a minimisation of least-squared error and visual fit to data is carried out based on

numerous model runs varying parameters that were identified in the sensitivity test.

# 5 Results and discussion

## 5.1 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⁻¹ (year 1), followed by an exponential growth phase to 46.8 µg C g⁻¹ (year 15), and then a decline to near-steady-state around 31.0 µg C g⁻¹, varying seasonally by roughly 12.0 µg C g⁻¹.
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Model output of the Damma Glacier forefield for a 75 year time series with nominal parameter values (Table 4).</figcaption></figure>

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_{\text{Sub}}$, $v_{\text{DIN}}$ and $v_{\text{DIP}}$ are adjusted. As with case study 1, optimisation is carried out based on the results of the sensitivity study, and a minimisation of least-squared error and visual fit to data is carried out based on

numerous model runs varying parameters that were identified in the sensitivity test.

# 5 Results and discussion

## 5.1 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⁻¹ (year 1), followed by an exponential growth phase to 46.8 µg C g⁻¹ (year 15), and then a decline to near-steady-state around 31.0 µg C g⁻¹, varying seasonally by roughly 12.0 µg C g⁻¹.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Sensitivity of model outputs ($\lambda$) to individual parameter variation. The model is forced with meteorological data from the Damma Glacier (Fig. 4) over 75 years.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Sensitivity of model outputs ($\lambda$) to individual parameter variation. The model is forced with meteorological data from the Damma Glacier (Fig. 4) over 75 years.</figcaption></figure>
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 \mu g \text{ C g}^{-1} \text{ yr}^{-1}$ (year 1) to $114.2 \mu g \text{ C g}^{-1} \text{ 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 \mu g \text{ C g}^{-1}$ (year 1) to $4129.2 \mu g \text{ 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.

## 5.2 Sensitivity and linearity

### 5.2.1 Sensitivity

Sensitivity analysis is presented in Fig. 8. The accumulation of autotrophic and heterotrophic biomass is most sensitive to variation in $Q_{10}$ ($\lambda \ge 0.70$). Biomass accumulation is also highly sensitive to adjustments in the active fraction ($d$) ($-0.55 \le \lambda \le -0.52$), the bioavailability of refractory substrate ($J_{S2}$) ($0.44 \le \lambda \le 0.64$), the partitioning of microbial necromass into labile and refractory pools ($q$) ($0.31 \le \lambda \le 0.45$) and microbial growth rates ($I_{\text{max}}$) ($-0.55 \le \lambda \le 0.67$). Biomass accumulation is moderately sensitive to death rates ($\alpha_A$ and $\alpha_H$), the efficiency of heterotrophic growth ($Y_H$) and the allochthonous substrate input ($v_{\text{Sub}}$) ($\pm 0.15 \le \lambda \le \pm 0.41$). Biomass accumulation has relatively low sensitivity ($\lambda \le \pm 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_{\text{Sub}}$ and $k_{\text{Sub}}$) and nitrogen fixers ($A_3$ and $H_3$) ($n_f, K_{N2}$), the bioavailability of labile substrate ($J_{S1}$), exudate rates ($\text{ex}_A$ and $\text{ex}_H$), the input of allochthonous nutrients ($v_{\text{DIN}}$ and $v_{\text{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 \le \lambda \le 0.00$ and $-0.05 \le \lambda \le 0.03$, respectively), but has a proportionally large effect on the accumulation of DIN, DIP and total nitrogen fixation ($0.22 \le \lambda \le 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 ($\lambda \le \pm 0.09$) but

strongly affects DIN, DIP and total nitrogen fixed ($\lambda = 0.96, \lambda = -0.96$ and $\lambda = 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_{\text{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_{\text{maxH}}$) leads to a substantial increase in both heterotrophic and autotrophic biomass ($\lambda = 0.47$ and $\lambda = 0.67$, respectively). However, strikingly, an increase in the maximum growth rate of autotrophs ($I_{\text{maxA}}$) has the net effect of lowering autotrophic and heterotrophic biomass ($\lambda = -0.36$ and $\lambda = -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_{\text{maxA}}$), their rapid growth increases nutrient scarcity ($\lambda_{\text{DIN}} = -0.57, \lambda_{\text{DIP}} = -1.58$) and heterotrophic growth becomes nutrient limited. However, if heterotrophs are given a competitive advantage (i.e. increasing $I_{\text{maxH}}$), substrate is degraded more rapidly, liberating DIN ($\lambda = 0.76$) and DIP ($\lambda = 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 ($\lambda = 0.47$)) in this oligotrophic and relatively nutrient poor environment. Similarly, increasing autotrophic death rates ($\alpha_A$) reduces competition for nutrient resources whilst also increasing the availability of degradable organic matter (autotrophic necromass), thereby increasing heterotrophic biomass ($\lambda = 0.24$).

An increase in the heterotrophic growth efficiency ($Y_H$) leads to an overall more rapid accumulation of heterotrophic biomass ($\lambda = 0.14$) at the expense of autotrophs ($\lambda = -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).
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 \mu g \text{ C g}^{-1} \text{ yr}^{-1}$ (year 1) to $114.2 \mu g \text{ C g}^{-1} \text{ 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 \mu g \text{ C g}^{-1}$ (year 1) to $4129.2 \mu g \text{ 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.

## 5.2 Sensitivity and linearity

### 5.2.1 Sensitivity

Sensitivity analysis is presented in Fig. 8. The accumulation of autotrophic and heterotrophic biomass is most sensitive to variation in $Q_{10}$ ($\lambda \ge 0.70$). Biomass accumulation is also highly sensitive to adjustments in the active fraction ($d$) ($-0.55 \le \lambda \le -0.52$), the bioavailability of refractory substrate ($J_{S2}$) ($0.44 \le \lambda \le 0.64$), the partitioning of microbial necromass into labile and refractory pools ($q$) ($0.31 \le \lambda \le 0.45$) and microbial growth rates ($I_{\max}$) ($-0.55 \le \lambda \le 0.67$). Biomass accumulation is moderately sensitive to death rates ($\alpha_A$ and $\alpha_H$), the efficiency of heterotrophic growth ($Y_H$) and the allochthonous substrate input ($v_{\text{Sub}}$) ($\pm 0.15 \le \lambda \le \pm 0.41$). Biomass accumulation has relatively low sensitivity ($\lambda \le \pm 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_{\text{Sub}}$ and $k_{\text{Sub}}$) and nitrogen fixers ($A_3$ and $H_3$) ($n_f$, $K_{N2}$), the bioavailability of labile substrate ($J_{S1}$), exudate rates ($\text{ex}_A$ and $\text{ex}_H$), the input of allochthonous nutrients ($v_{\text{DIN}}$ and $v_{\text{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 \le \lambda \le 0.00$ and $-0.05 \le \lambda \le 0.03$, respectively), but has a proportionally large effect on the accumulation of DIN, DIP and total nitrogen fixation ($0.22 \le \lambda \le 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 ($\lambda \le \pm 0.09$) but

strongly affects DIN, DIP and total nitrogen fixed ($\lambda = 0.96$, $\lambda = -0.96$ and $\lambda = 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_{\text{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_{\max H}$) leads to a substantial increase in both heterotrophic and autotrophic biomass ($\lambda = 0.47$ and $\lambda = 0.67$, respectively). However, strikingly, an increase in the maximum growth rate of autotrophs ($I_{\max A}$) has the net effect of lowering autotrophic and heterotrophic biomass ($\lambda = -0.36$ and $\lambda = -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_{\max A}$), their rapid growth increases nutrient scarcity ($\lambda_{\text{DIN}} = -0.57$, $\lambda_{\text{DIP}} = -1.58$) and heterotrophic growth becomes nutrient limited. However, if heterotrophs are given a competitive advantage (i.e. increasing $I_{\max H}$), substrate is degraded more rapidly, liberating DIN ($\lambda = 0.76$) and DIP ($\lambda = 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 ($\lambda = 0.47$)) in this oligotrophic and relatively nutrient poor environment. Similarly, increasing autotrophic death rates ($\alpha_A$) reduces competition for nutrient resources whilst also increasing the availability of degradable organic matter (autotrophic necromass), thereby increasing heterotrophic biomass ($\lambda = 0.24$).

An increase in the heterotrophic growth efficiency ($Y_H$) leads to an overall more rapid accumulation of heterotrophic biomass ($\lambda = 0.14$) at the expense of autotrophs ($\lambda = -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).
## 5.2.2 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$, $\alpha_A$, $\alpha_H$, $v_{DIN}$ and $v_{DIP}$. Therefore, the sensitivity value $\lambda$ is likely to change when looking at values away from the nominal parameter value, as illustrated for $I_{maxA}$ in Fig. 6b. For non-linear parameters $\alpha_A$, $\alpha_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.

## 5.3 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 $\alpha_A$, $\alpha_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.,

<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

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
## 5.2.2 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$, $\alpha_A$, $\alpha_H$, $v_{DIN}$ and $v_{DIP}$. Therefore, the sensitivity value $\lambda$ is likely to change when looking at values away from the nominal parameter value, as illustrated for $I_{maxA}$ in Fig. 6b. For non-linear parameters $\alpha_A$, $\alpha_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.

## 5.3 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 $\alpha_A$, $\alpha_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.,

<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

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
<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

into models because of the assumptions and simplifications in model design.

Microbial death rates ($\alpha_A$ and $\alpha_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 ($\lambda \le -0.13$ and $\lambda \ge 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 \text{ 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_{\text{Sub}}$, $v_{\text{DIN}}$ and $v_{\text{DIP}}$ are wide, and the resulting uncertainty is very high (75 to 2503 %). Furthermore, variation in parameters $v_{\text{DIN}}$ and $v_{\text{DIP}}$ results in highly non-linear behaviour, with comparatively large changes resulting from small changes in the parameter value. However, the effect of $v_{\text{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 \le \lambda \le 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 \mu\text{g N g}^{-1} \le K_{N2} \le 1128 \mu\text{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}$ ($\lambda = 0.08$ and $\lambda = 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 ($\alpha_A$ and $\alpha_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 ($\text{ex}_A$ and $\text{ex}_H$) (uncertainty $\le 2\%$) have relatively tight parameter bounds (Allison, 2005), as well as low overall sensitivity ($-0.01 \le \lambda \le$
<figure><img src="image_15.png" /><figcaption>**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.</figcaption></figure>

into models because of the assumptions and simplifications in model design.

Microbial death rates ($\alpha_A$ and $\alpha_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 ($\lambda \le -0.13$ and $\lambda \ge 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 \text{ 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_{\text{Sub}}$, $v_{\text{DIN}}$ and $v_{\text{DIP}}$ are wide, and the resulting uncertainty is very high (75 to 2503 %). Furthermore, variation in parameters $v_{\text{DIN}}$ and $v_{\text{DIP}}$ results in highly non-linear behaviour, with comparatively large changes resulting from small changes in the parameter value. However, the effect of $v_{\text{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 \le \lambda \le 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 \mu\text{g N g}^{-1} \le K_{N2} \le 1128 \mu\text{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}$ ($\lambda = 0.08$ and $\lambda = 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 ($\alpha_A$ and $\alpha_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 ($\text{ex}_A$ and $\text{ex}_H$) (uncertainty $\le 2\%$) have relatively tight parameter bounds (Allison, 2005), as well as low overall sensitivity ($-0.01 \le \lambda \le$
<figure><img src="image_15.png" /><figcaption>Figure 11. Model output optimised to observational data from the Damma Glacier, Switzerland.</figcaption></figure>

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 $\mu 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 ($\lambda = -0.07$), as is reflected in other models (Ingwersen et al., 2008; Blagodatsky and Richter, 1998).

## 5.4 Optimisation and model dynamics

### 5.4.1 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 $\mu 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 $\mu g \ C \ g^{-1}$ (years 30 to 75). The microbial community

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. Labile substrate is
<figure><img src="image_16.png" /><figcaption>Figure 11. Model output optimised to observational data from the Damma Glacier, Switzerland.</figcaption></figure>

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 $\mu g \text{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 ($\lambda = -0.07$), as is reflected in other models (Ingwersen et al., 2008; Blagodatsky and Richter, 1998).

## 5.4 Optimisation and model dynamics

### 5.4.1 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 \mu g \text{C g}^{-1}$ (year 12), followed by a decline over the following 10 years, after which biomass is fairly stable at roughly $37.0 \mu g \text{C g}^{-1}$ (years 30 to 75). The microbial community

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. Labile substrate is
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Evolution of the microbial community composition in model output optimised from the Damma Glacier, Switzerland, over the first 20 years of soil formation.</figcaption></figure>

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⁻¹ yr⁻¹) 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⁻¹ yr⁻¹) 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⁻¹ yr⁻¹ (Schulz et al., 2013; Guelland et al., 2013a), which is within the range predicted by the SHIMMER 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⁻¹ yr⁻¹), carbon substrate accumulates to roughly half of the highest maximum substrate sampled in the BigLink project (31 363.1 µg C g⁻¹) (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 (Antarcitica) (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⁻¹) 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⁻¹ yr⁻¹ (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
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Evolution of the microbial community composition in model output optimised from the Damma Glacier, Switzerland, over the first 20 years of soil formation.</figcaption></figure>

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⁻¹ yr⁻¹) 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⁻¹ yr⁻¹) 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⁻¹ yr⁻¹ (Schulz et al., 2013; Guelland et al., 2013a), which is within the range predicted by the SHIMMER 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⁻¹ yr⁻¹), carbon substrate accumulates to roughly half of the highest maximum substrate sampled in the BigLink project (31 363.1 µg C g⁻¹) (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 (Antarcitica) (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⁻¹) 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⁻¹ yr⁻¹ (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
<figure><img src="image_17.png" /><figcaption>**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.</figcaption></figure>

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 \mu g C g^{-1} d^{-1}$) than winter ($0.05 \mu 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).

### 5.4.2 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 \mu g C g^{-1}$ in year 1 to $11.5 \mu g C g^{-1}$ in year 40, and at roughly $0.02 \mu 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_{1-3}$) make up 17.2% of the biomass in year 15. Soil heterotrophs ($H_2$) account for only 4.7% of biomass in year 15 ($0.33 \mu g C g^{-1}$), and are severely limited by DIN
<figure><img src="image_18.png" /><figcaption>**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.</figcaption></figure>

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 \mu g C g^{-1} d^{-1}$) than winter ($0.05 \mu 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).

### 5.4.2 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 \mu g C g^{-1}$ in year 1 to $11.5 \mu g C g^{-1}$ in year 40, and at roughly $0.02 \mu 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_{1-3}$) make up 17.2% of the biomass in year 15. Soil heterotrophs ($H_2$) account for only 4.7% of biomass in year 15 ($0.33 \mu g C g^{-1}$), and are severely limited by DIN
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Model output optimised to observational data from the Athabasca Glacier, Canada.</figcaption></figure>

(11.3 µg N g⁻¹); 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.

## 6 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.

### 6.1 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). Resolution of autotrophic and heterotrophic biomass would help in validating model predictions of community