A parameterization of respiration in frozen soils based on substrate availability

Respiration in frozen soils is limited to thawed substrate within the thin water films surrounding soil particles. As temperatures decrease and the films become thinner, the available substrate also decreases, with respiration effectively ceasing at −8 C. Traditional exponential scaling factors to model this effect do not account for substrate availability and do not work at the century to millennial timescales required to model the fate of the nearly 1100 Gt of carbon in permafrost regions. The exponential scaling factor produces a false, continuous loss of simulated permafrost carbon in the 20th century and biases in estimates of potential emissions as permafrost thaws in the future. Here we describe a new frozen biogeochemistry parameterization that separates the simulated carbon into frozen and thawed pools to represent the effects of substrate availability. We parameterized the liquid water fraction as a function of temperature based on observations and use this to transfer carbon between frozen pools and thawed carbon in the thin water films. The simulated volumetric water content (VWC) as a function of temperature is consistent with observed values and the simulated respiration fluxes as a function of temperature are consistent with results from incubation experiments. The amount of organic matter was the single largest influence on simulated VWC and respiration fluxes. Future versions of the parameterization should account for additional, non-linear effects of substrate diffusion in thin water films on simulated respiration. Controlling respiration in frozen soils based on substrate availability allows us to maintain a realistic permafrost carbon pool by eliminating the continuous loss caused by the original exponential scaling factors. The frozen biogeochemistry parameterization is a useful way to represent the effects of substrate availability on soil respiration in model applications that focus on century to millennial timescales in permafrost regions.


Introduction
Incubated frozen soil samples show a strong decrease in respiration with temperature below freezing associated with decreased substrate availability (Mikan et al., 2002).Here we define substrate availability as the portion of soil organic matter accessible as a medium for microbial activity.The sharp decline in respiration results from the fact that soils maintain some liquid water at temperatures below freezing (Romanovsky and Osterkamp, 2000;Davis, 2001;Kurylyk and Watanabe, 2013).The mechanism is well understood: freezing point depression driven by the curvature of water around small, hydrophilic soil particles, analogous to freezing point depression caused by solutes in water (Davis, 2001).The result is thin liquid water films surrounding soil particles at temperatures below freezing.Essentially, microbial activity becomes limited to available thawed organic matter or dissolved organic carbon (DOC) within the thin water films.The diffusion of substrates to the microbes becomes limited to narrow water channels in the thin water films (Rivkina et al., 2000).As temperatures drop below freezing, the water films become thinner and the available substrate and associated microbial activity rapidly decreases, with respiration effectively ceasing below temperatures of −7 to −8 • C (Oechel et al., 1997;Mast et al., 1998;Hobbie et al., 2000;Mikan et al., 2002).
The most common way to model this sharp decline in respiration below 0 • C is to apply an exponential temperature scaling factor to the simulated microbial decay rates: where S f is a freezing scaling factor, T is soil temperature, T ref is a reference temperature (typically 0 • C), and Q 10f is the change in respiration for a 10 K drop in temperature below freezing.Substrate availability in frozen soil is determined by the amount of thawed organic matter and the diffusion of DOC in the thin water films.The Q 10f formulation does not account for substrate availability, but rather is based on the Arrhenius equation for kinetic controls on respiration in thawed soils: where S T is a thawed respiration scaling factor, T ref is a reference temperature (typically 5 or 10 • C), and Q 10 is the change in respiration changes for a 10 • C change in temperature (Raich and Schlesinger, 1992;Denning et al., 1996;Potter et al., 1993).The kinetic Q 10 formulation is based on the Arrhenius equation for chemical reaction rates: the higher the temperature, the greater the number of molecules that have energies greater than the minimum activation energy and the faster the rate of microbial decay.Q 10 varies between 1.5 and 3.0 based on incubation studies or analysis of eddy covariance flux data (Oechel et al., 1997;Mast et al., 1998;Hobbie et al., 2000;Mikan et al., 2002).In contrast, Q 10f varies between 164 and 237 based on incubation of frozen soil samples (Mikan et al., 2002).
Because substrate availability rather than kinetics controls respiration below freezing, the Q 10f formulation is inappropriate when attempting to model the large amount of frozen organic matter in permafrost regions.Permafrost is soil at or below 0 • C for at least 2 years and permafrost soils in the high northern latitudes contain ∼ 1100 Gt of carbon, mostly in the form of frozen organic matter (Tarnocai et al., 2009;Hugelius et al., 2013).This frozen permafrost carbon was buried over millennial timescales by alluvial sedimentation, dust deposition, and peat development, which increase soil depth and freeze organic matter at the bottom of the active layer into the permafrost (Zimov et al., 2006a, b;Schuur et al., 2008;Ping et al., 2015).Permafrost carbon was buried during or since the last ice age with maximum ages between 15 000 and 30 000 years (Schuur et al., 2008;Dutta et al., 2006).Most models do not include these burial mechanisms and simply initialize the frozen carbon based on observed carbon densities in permafrost (Schaefer et al., 2011).Applying the Q 10f formulation results in effective turnover times of 200 to 500 years such that the simulated permafrost carbon does not persist long enough in the model to match observed vertical carbon distributions (Harden et al., 2012;Hugelius et al., 2013).To counter this slow loss of permafrost carbon, Koven et al. (2011) increased Q 10f to 1000 and we initially increased Q 10f to 10 000, which are well beyond observed values.These large Q 10f values increased the effective turnover time in permafrost to 10 000 years or more and allowed the proper buildup of permafrost carbon.However, they also had the undesired effect of effectively shutting down wintertime respiration, which can account for a significant fraction of total annual respiration (Alm et al., 1999;Fang et al., 1999;Zimov et al., 1993aZimov et al., , b, 1996;;Schmidt and Lipson, 2004).The problem is that the Q 10f formulation based on the Arrhenius equation does not account for substrate availability, making it inappropriate when representing permafrost carbon dynamics on timescales of 500 to 10 000 years.
Here we present a new frozen biogeochemistry parameterization based on substrate availability.We link the physical processes that determine the amount of liquid water in frozen soils to a fully prognostic set of biogeochemical carbon pools.Tucker (2014) successfully combined liquid water content in frozen soils with an enzyme kinetic model of respiration accounting for DOC diffusion.We integrated the liquid water content of frozen soils into the biogeochemistry parameterization of the Simple Biosphere/Carnegie-Ames-Stanford Approach (SiBCASA) model.The frozen biogeochemistry parameterization separates kinetic controls from substrate availability in frozen soils to simultaneously reproduce observed liquid water fractions and frozen carbon content in permafrost.

Methods
SiBCASA is a full land surface parameterization with fully integrated water, energy, and carbon cycles (Schaefer et al., 2008).SiBCASA predicts the moisture content and temperature of the canopy, canopy air space, and soil (Sellers et al., 1986;Vidale and Stockli, 2005).Schaefer et al. (2009) modified the snow module to better simulate permafrost dynamics and SiBCASA has been used to predict future permafrost degradation and global carbon emissions from thawing permafrost (Schaefer et al., 2011).The SiBCASA soil model has 25 layers to a depth of 15 m with geometrically increasing layer thickness with depth starting with 2 cm at the surface.SiBCASA simultaneously predicts soil moisture and temperature for each soil layer on a 10 min time step.The soil moisture model accounts for surface infiltration, surface runoff, and subsurface runoff.The soil model separately tracks liquid and frozen water at each time step as prognostic variables, accounting for the effects of latent heat (Schaefer et al., 2008(Schaefer et al., , 2009)).The soil thermodynamic and hydraulic properties are a volume-weighted fraction of the properties of liquid water, ice, mineral soil, and organic matter (Schaefer et al., 2009).SiBCASA does not include any gas diffusion processes within soil pore spaces and does not include solute diffusion processes in the soil water.
SiBCASA predicts soil organic matter, surface litter, and live biomass (leaves, roots, and wood) in a system of 13 prognostic carbon pools (Schaefer et al., 2008).This includes four live pools (starch, leaves, roots, and wood) and nine soil carbon pools (coarse woody debris, slow, etc.).SiBCASA does not currently include a DOC pool.SiBCASA represents the biogeochemistry as a system of coupled, first-order linear differential equations: where C i is the ith carbon pool, t is time, S T and S f are the temperature and freezing scaling factors as described above, S W is a soil moisture scaling factor, τ is the reference turnover time for the pool, and G represents gains from other carbon pools (Schaefer et al., 2008).The first term represents the decay of organic material, some fraction of which is released as respiration and the rest transferred to other pools, from the coarse woody debris to the slow pool, for example (see Schaefer et al., 2008, for a complete description).
SiBCASA uses a Q 10 of 1.5, which is held constant for all temperatures.Prior to incorporating the frozen biogeochemistry parameterization, SiBCASA used a Q 10f of 200, which was also held constant, but only for temperatures below 0 • C. Each soil layer has a complete set of prognostic soil carbon pools.Once per day for each soil layer SiBCASA recalculates the organic soil fraction (f org ) used to determine thermodynamic and hydraulic properties.SiBCASA redistributes carbon in the top soil layers to simulate the surface organic layer (Jafarov and Schaefer, 2016).This allows for the buildup of a surface organic layer crucial to simulating soil thermodynamics in permafrost regions.However, SiB-CASA does not represent the cryoturbation and sedimentation processes required to build up a large permafrost carbon pool.Instead, we initialize the permafrost carbon content using the observed distribution from the Northern Circumpolar Soil Carbon Database version 2 (NCSCDv2) (Hugeluis et al., 2013).
Substrate availability is determined by the portion of organic matter that is thawed and by the diffusion of DOC in the thin water films.SiBCASA does not have a DOC pool and does not include DOC diffusion, so the new frozen biogeochemistry parameterization represents only the effects of the amount of thawed organic matter on substrate availability.The frozen biogeochemistry parameterization uses three sets of carbon pools for each soil layer: thawed, frozen film, and bulk frozen pools (Fig. 1).The bulk frozen pools represent frozen carbon in soil pore spaces away from the thin liquid water films surrounding mineral soil particles.The frozen film pools represent frozen soil carbon within the maximum extent of the thin water films.The thawed pools represent organic matter in the thin liquid water films.Soil carbon in the bulk and film frozen pools is unavailable for microbial decay and remain static and sequestered until thawed.Simulated microbial decay and associated respiration occurs only A schematic of the new frozen biogeochemistry parameterization around a single, idealized soil particle.Soil carbon is divided into thawed carbon in the thin water films surrounding soil particles, film frozen carbon in the maximum extent of the thin liquid water film, and bulk frozen carbon in the pore spaces between soil particles.The thawed carbon represents available substrate for metabolism by microbes in the thin water film.
in the thawed carbon pools in the thin water films and will continue to decay below 0 • C.This complete separation of frozen and thawed soil carbon represents the effect of the amount of thawed organic matter on substrate availability in the thin liquid water films for microbial metabolism and respiration in frozen soils.This new frozen carbon parameterization replaced the original Q 10f formulation.
The thawed, film frozen, and bulk frozen carbon pools all have the same 13-pool structure as defined in Schaefer et al. (2008).The frozen biogeochemistry parameterization applies only to the nine soil carbon pools.The live biomass pools (starch, leaves, fine roots, and wood) are always considered "thawed" because the growth and mortality processes that govern them do not depend upon soil texture and associated thin water films.Carbon is transferred between the frozen and corresponding thawed pools as the amount of liquid water changes with simulated soil temperature: thawed "slow" pool to bulk frozen "slow" pool, etc.
We assumed that the mineral soil, carbon, liquid water, and ice are uniformly distributed within the soil layer such that the liquid water fraction equals the thawed fraction.We calculate the liquid water fraction for each soil type using a modified power law formulation introduced by Lovell (1957) and refined by Nicolsky et al. (2009): where ϕ i is the fraction of liquid water for the ith soil type, T is the soil layer temperature (  ϕ i ranging from simple piecewise continuous linear models to full physics models based on the Clapeyron equation.We desired a formulation that produced realistic results, but was not overly complicated, so we decided upon the power law formulation.We assumed four soil types, pure clay, silt, sand, and organic matter, each with different values of b i (Table 1).
We used b i values from Romanovsky and Osterkamp (2000) and Kurylyk and Watanabe (2013).We calculated the ϕ crit by inserting 0 • C into Eq.( 4) above.Clay represents the finest grained soil material with the greatest amount of liquid water below freezing and organic material represents the coarsest soil material with the least amount of liquid water.
We made two important changes to the Nicolsky et al. (2009) model in order to reproduce the observed ϕ i discontinuity at 0 • C and differentiate between the bulk frozen and film frozen pools.The original Nicolsky et al. (2009) formulation produced a continuous ϕ i for all values of T , but liquid water observations show a discontinuity at 0 • C where the actual ϕ i is determined by bulk liquid to ice conversion (Davis, 2001).First, SiBCASA does not include aqueous chemistry to calculate the freezing point depression due to dissolved solutes, so we changed T * from freezing point depression to a simple offset temperature.Second, we fixed T ref at 0.1 • C rather than allowing it to change with depth and temperature (Nicolsky et al., 2009).Together, these two changes reproduced the observed ϕ i discontinuity at 0 • C. ϕ crit represents the liquid water content at 0 • C and determines the boundary between bulk freezing and thin film processes.Bulk freezing dominated by latent heat effects occurs for ϕ>ϕ crit , while thin film effects dominate for ϕ ≤ ϕ crit .Essentially, ϕ crit represents the maximum thickness of the thin water films around soil grains and allows us to differentiate between the bulk frozen and film frozen pools.
Figure 2 shows the parameterized ϕ i as a function of temperature for clay, silt, sand, and organic soils.The thin dashed lines represent ϕ crit for each soil type defining the boundary between bulk freezing and thin film effects.The shapes and magnitudes of the individual curves are consistent with observed ϕ values (Davis, 2001).The liquid water fraction below freezing varies between less than 1 % to as high as 50 %, depending on soil texture: larger soil particles have more negative values of b i and less liquid water when frozen.For organic soil, we assumed the individual particles of organic matter were, on average, too large to produce unfrozen water and assigned the largest b i resulting in the smallest val- ues of ϕ i .ϕ i is next lowest in frozen sandy soils because the particles are large and hydrophobic with little unfrozen water.ϕ i is highest in frozen clay soils because the particles are small and hydrophilic with the largest amount of unfrozen water.
We assumed that the overall liquid water fraction at any given temperature is the weighted average of those for pure clay, silt, sand, and organic matter: where ϕ is the liquid water fraction for each soil layer, f org is the volumetric organic soil fraction, and f c , f si , and f sa are the volumetric fractions of clay, silt, and sand for the mineral portion of the soil.ϕ c , ϕ si , ϕ sa , and ϕ o are the liquid water fractions for pure clay, silt, sand, and organic matter as a function of T based on the power law formulation above.f org is calculated as the ratio of simulated carbon density to the observed density of pure organic matter: where M o is the simulated organic matter mass per soil layer, z is the thickness of the soil layer, and ρ o is the observed bulk density of pure organic soil, which we assumed was 140 kg m −3 (Schaefer et al., 2009).M o is the sum of all prognostic carbon pools in each soil layer, converted from mass of pure carbon to mass of organic matter assuming organic matter is 50 % carbon (Schaefer et al., 2009).The calculation of ϕ crit is also a weighted average of the ϕ crit values for each soil type.The carbon content of each soil layer varies relatively slowly with time depending on the prognostic carbon pools and we assume that the soil matrix and associated physical properties do not change when frozen.Consequently, SiB-CASA recalculates f o and ϕ crit for each soil layer once per day for T ≥ 0 • C. The clay, silt, and sand fractions of the mineral portion of the soil are constant in time, but vary in Figure 3.A schematic showing the transfers of carbon between the thawed, bulk frozen, and film frozen as the liquid water fraction (ϕ) changes with time and temperature (T ).As the soil freezes, ϕ follows the red line and carbon is transferred from the thawed pool to the bulk frozen carbon pool until ϕ reaches ϕ crit , then carbon is transferred from the thawed to the film frozen pool.The reverse happens during thawing.
space based on the Harmonized World Soil Database (Wei et al., 2014).
SiBCASA transfers carbon between the thawed, film, and bulk frozen pools each time step depending on the change in ϕ: φ = φ t − φ t−1 , where ϕ is the change in thawed soil fraction for a single time step, the superscript t denotes the value at the current time step, and t − 1 the value at the previous time step (Fig. 3).A negative ϕ indicates freezing soil and a positive ϕ indicates thawing soil.ϕ stays constant at 1.0 until the soil layer reaches 0 • C. As freezing begins, T stays constant at 0 • C and ϕ decreases to account for the latent heat of fusion for water.When ϕ reaches ϕ crit , T decreases below 0 • C and ϕ follows the liquid water curve.During thaw, ϕ follows the reverse path, ignoring possible hysteresis effects.ϕ>ϕ crit represents freezing or thawing of the bulk pore space between soil grains and ϕ ≤ ϕ crit represents freezing and thawing of the thin films around soil grains.During freezing, the bulk carbon away from the soil grains freezes first and the liquid water film freezes last.During thaw, the reverse is true with the frozen film carbon thawing first and the bulk carbon last.
Table 2 shows the carbon bookkeeping rules governing the transfer of carbon between thawed, film frozen, and bulk frozen carbon pools each time step.C i_thaw is the ith thawed carbon pool, C i_film is the ith film frozen carbon pool, and C i_bulk is the ith bulk frozen pool, where i represents the nine soil carbon pools described above.δ i_2 bulk is the carbon transferred from the ith thawed carbon pool to the ith bulk frozen pool, with similar nomenclature for transfers to and from the film frozen pools.Table 2 is organized by freezing starting with bulk pools first and film frozen pools second, with the reverse for thawing.When ϕ crosses ϕ crit , part of the carbon goes to the bulk pools and part to the film frozen pools.The same equations are applied in sequence to all nine soil carbon pools each time step.
We compared simulated response curves for volumetric water content (VWC) vs. temperature against observed values at Bonanza Creek, Alaska (Jafarov et al., 2013).Simulated VWC is ϕθ η, where θ is the soil saturation fraction and η is soil porosity for a given soil layer.θ is the fraction of pore space filled with both liquid water and ice.η varies from 0.85 for pure organic matter to between 0.35 and 0.45 for pure mineral soil, depending on soil texture.We used VWC data collected at the Bonanza Creek Long Term Ecological Research site in central Alaska using a Hydro Probe soil moisture sensor from 2009 to 2014.Bonanza Creek is in the discontinuous permafrost zone and we used VWC data from an unburned and recently burned site.Both sites have sensors installed at 19, 36, and 54 cm depths (Romanovsky and Osterkamp, 2001;Romanovsky et al., 2003).At the unburned site, all three depths fall in the surface organic layer, while at the burned site, all three depths fall within mineral soil below the organic layer.Because of the huge influence of f org on ϕ and thus VWC, we separately compared simulated to observed values for organic soil and mineral soil.The simulated organic layer thickness was 21 cm, so we extracted the simulated VWC in the organic layer at 19 cm depth and compared them with organic soil observations at 19 cm depth at the unburned site.We extracted simulated VWC for mineral soil at 54 cm depth and compared them with mineral soil observations at 54 cm depth at the burned site.
We compared response curves for simulated respiration vs. temperature against observed values from four independent incubations of frozen soil samples: Rivkina et al. (2000), Mikan et al. (2002), Eberling and Brandt (2003), and Panikov et al. (2006).Each of these studies collected samples from different locations and used different temperature ranges, durations, and protocols for the incubations (Table 3).We converted the observed respiration values to common units of flux per mass of carbon per day (µg C g −1 C d −1 ) using total organic carbon (TOC) values and soil bulk densities as noted in each paper.TOC is the ratio of organic matter mass to total dry soil mass, which we converted to f org to help evaluate model output.Rivkina et al. (2000) collected 14 C counts per minute, which we could not convert into a respiration flux, so we only did a qualitative comparison.Panikov et al. (2006) did not report an observed TOC, so we estimated the TOC from the observed bulk densities.For our various unit conversions, we assumed a mineral soil density of 1400 kg m −3 and a ρ o of 140 kg m −3 (Schaefer et al., 2009).We converted the normalized values from Eberling and Brandt (2003) to absolute values using the observed respiration at the reference temperature of 7 • C. We averaged the Mikan et al. (2002) results for individual soil samples to get an average respiration at each incubation temperature, consistent with reported values in Elberling and Brandt (2003) and Panikov et al. (2006).
We compared results for the frozen biogeochemistry parameterization using two point simulations: Toolik Lake on www.biogeosciences.net/13/1991/2016/Biogeosciences, 13, 1991-2001, 2016 Thawing ϕ>0 ϕ t−1 <ϕ crit ϕ t >ϕ crit film and bulk to thaw We compared the Toolik Lake simulation output against the incubation data and the Bonanza Creek output against the VWC data.We ran simulations at Bonanza Creek and all the sites in Table 3, but because we are evaluating the simulated temperature response functions of VWC and respiration, the comparisons at the various sites became extremely repetitive.
Here we include only the Toolik Lake and Bonanza Creek results because they cover a significant portion of the temperature ranges in the data.The location of the simulation had almost no effect on our results, as long as the simulated soil temperature ranges overlapped sufficiently with the incubation temperatures.More than any other factor, the simulated f org dominated the simulated VWC and respiration response functions.So rather than a site-by-site comparison, we focused on a comparison of organic and inorganic soils by choosing SiBCASA soil layers either within the simulated surface organic layer or below it.The results for the remaining sites are nearly identical to the Toolik Lake and Bonanza Creek simulations.
As input weather, we used the Climatic Research Unit National Center for Environmental Predictions (CRUNCEP) data set (Wei et al., 2014), extracting the CRUNCEP data for the grid cell containing each point.CRUNCEP is reanalyzed weather data at 0.5 × 0.5 degree latitude and longitude resolution optimally consistent with a broad array of observations spanning 110 years, from 1901 to 2010.Starting from steady-state initial conditions, we ran the point simulations from 1901 to 2010 and extracted model output for the 5 years closest to the time period covered by the observations.The Bonanza Creek observations cover 2009 to 2014, so we extracted model output from 2005 to 2010.This slight mismatch is reasonable since our objective was to evaluate the simulated temperature response functions rather than compare simulated and observed time series.Mikan et al. (2002) collected the soil samples in 1998, so we extracted model output from 1996 to 2000 for the Toolik Lake simulation.
We also ran simulations with and without the new frozen biogeochemistry parameterization for the entire permafrost domain in the Northern Hemisphere (Brown et al., 1997).We initialized the permafrost carbon content within the top 3 m of permafrost using the observed distribution from the Northern Circumpolar Soil Carbon Database version 2 (NC-SCDv2) (Hugeluis et al., 2013).The permafrost carbon was split between the bulk frozen and film frozen pools based on the ϕ and ϕ crit values at the start of spinup simulation.We selected the first 30 years from the CRUNCEP data set (1901 to 1931) and randomly distributed them over 4000 years to spin SiBCASA up to steady-state initial conditions.We spun the model up to steady-state initial conditions in 1900 with the new frozen biogeochemistry parameterization turned on.We then ran two simulations from 1901 to 2010 starting from the same initial conditions, one with the new frozen biogeochemistry parameterization and one without.

Results
Our ϕ parameterization produced a VWC vs. temperature response consistent with observed values (Fig. 4).Below freezing, the simulated VWC for organic soil is higher than observed values while the simulated VWC for mineral soil closely matches observed values.Above freezing, both the observed and simulated VWC for organic soil varied widely because of a strong variation in simulated saturation fraction over time.SiBCASA assumes a porosity of 0.85 for organic soil, which results in a VWC at saturation just below the maximum observed values of ∼ 0.9.For mineral soil, the simulated and observed VWC both stay near saturation, but the simulated values are greater than observed above freezing.This difference above freezing probably results from the assumed soil texture in our simulation, but we lacked observations of soil texture at Bonanza Creek to confirm this.
The frozen biogeochemistry parameterization reproduced the VWC discontinuity at 0 • C, but the observed VWC shows some noise at 0 • C because we converted to daily averages.A large observed diurnal cycle in soil temperature resulted in conditions where the soil froze at night and thawed during the day.Thus, the daily average temperature may be 2 • C, for example, but the daily average VWC reflects both frozen and unfrozen conditions.This produces a horizontal "spread" in the VWC around 0 • C, and thus noise in the discontinuity.The SiBCASA soil model had less diurnal variability in simulated soil temperature, so the noise was less than the observed values.The spread appears in observed values for both the mineral and organic soils, but is greater for the mineral soil because the observed temperature amplitude is greater.
Although each incubation study sampled from different locations and used different protocols, the overall results show consistent magnitudes and patterns (Fig. 5).The observed respiration values for Plotnikovo and Koppara were so close that we plotted them using the same color for clarity.Respiration decreased exponentially with decreasing temperature, with a faster rate of decline below freezing.The observed respiration at 0 • C is ∼ 60 µg C g −1 C d −1 , with no obvious discontinuity.The respiration decreases to ∼ 2 µg C g −1 C d −1 at −10 • C, a reduction of 97 %.Observed respiration below −10 • C varies between 0.1 and 2 µg C g −1 C d −1 representing a 97-99 % reduction compared to 0 • C, with the Barrow incubations showing residual respiration at temperatures as low as −39 • C. Except for Plotnikovo and Koppara, all the incubated soil samples had TOC <0.07 (f org <0.45), indicating a mix of mineral and organic soil with corresponding higher values of ϕ and respiration.In contrast, the Plotnikovo and Koppara samples were pure organic matter with TOC = 1.0 (f org = 1.0) and, because ϕ was much less, showed the lowest observed respiration of all the studies.
Above freezing, the simulated respiration as a function of temperature was consistent in magnitude with observed values from incubation experiments, but showed much greater variability (Fig. 6).The variability in simulated respiration  3. The Koppara and Plotnikovo incubation results appear as the same color.
above freezing results from temporal variability in simulated soil moisture content.Of course, the incubation data show no such variability because experiment protocols held soil moisture constant.Organic soil has greater hydraulic conductivity and porosity than mineral soil with corresponding greater variability in soil moisture, and thus respiration.On average, the simulated respiration for organic soil above freezing was greater than the incubation data because of higher values of simulated TOC than in the soil samples.For the mineral soil, the average simulated respiration above freezing is consistent with the incubation values.
At 0 • C, temporal variations in simulated hydraulic conductivity produced a "tail" of low simulated respiration values in the mineral soil simulation results.The effective pore space and associated hydraulic conductivity decreases as liquid water changes to ice during the freezing process.This produced an increase in VWC and an associated increase in respiration as the simulated soil froze during fall and early winter.The effect occurs for both organic and mineral soil, but is more prominent in mineral soil because it is deeper in the soil column and because thawed mineral soil has lower hydraulic conductivity than thawed organic soil.Again, the incubation data do not show such an effect because soil moisture is held constant.
Below freezing, the simulated organic soil has lower respiration than mineral soil (Fig. 6a), consistent with the Plotnikovo and Koppara observations and demonstrating the strong influence of organic matter on frozen soil respiration.The observed TOC for nearly all of the incubation samples varied between 0.01 and 0.07, consistent with the simulated TOC of 0.04 for the mineral soil (f org = 0.27).The simulated TOC for the organic soil is 0.7 (f org = 0.95), much closer to the observed TOC of 1.0 for the Plotnikovo and Koppara incubations.The high TOC results in lower ϕ and respiration such that the simulated respiration matches the Plotnikovo and Koppara incubation data within 10 % at all temperatures.In contrast, the simulated respiration from mineral soil was less than the incubation data between 0 and −5 • C, and greater than the incubation data for temperatures less than −8 • C (Fig. 6b).The frozen biogeochemistry parameterization produced a sharp discontinuity in simulated respiration at 0 • C that mirrors the VWC discontinuity.The incubation data do not indicate such a sharp discontinuity, resulting in lower respiration than observed for temperatures between 0 and −5 • C.
The frozen biogeochemistry parameterization maintains a realistic permafrost carbon pool in the model.Figure 7 shows the total simulated permafrost carbon in the top 3 m of soil for the entire Northern Hemisphere permafrost domain, both with and without the frozen biogeochemistry parameterization.Hugelius et al. (2014) indicate a total of 800 Gt of carbon frozen in permafrost, with 619 Gt in the top 3 m, or 10 % higher than the 560 Gt C we simulate.Using the old kinetic Q 10f formulation, the permafrost carbon decreased at a nearly linear rate as it slowly decayed.In contrast, the new frozen biogeochemistry parameterization based on substrate availability allows SiBCASA to maintain a nearly constant permafrost carbon pool throughout most of the 20th century.After 1950, the temperatures slowly rise, resulting in thawing of ∼ 3 Gt of permafrost carbon representing 0.6 % of the initial pool.

Discussion
We hypothesize that the diffusivity of DOC and microbial waste products, and thus respiration, does not respond linearly to changes in VWC.Because we use the ϕ curve to explicitly define substrate availability, the frozen biogeochemistry parameterization assumes a linear response to VWC.As a result, the shape of the simulated respiration curve exactly matches the simulated VWC curve, with a discontinuity in simulated respiration at 0 • C not seen in the incubation data.As the bulk of the water in the soil freezes, the DOC concentration in the thin water films increases.This counteracts the decrease in the amount of thawed organic matter and decline in DOC diffusivity as the thickness of the thin films decrease from 15 to 5 nm between −1.5 and −10 • C (Rivkina et al., 2000;Tucker, 2014).The result is a nonlinear response in respiration to changing VWC between 0 and −5 • C. Rivkina et al. (2000) found observed 14 C counts in respired CO 2 exactly matched the VWC curve, but they infused their samples with 14 C-labeled glucose, which may not be strongly affected by decreases in diffusivity.In contrast, Panikov et al. (2006) achieved an R 2 of 0.99 with an exponential curve fit of respiration to both temperature and VWC.Since VWC is a function of soil temperature, this additional dependency on VWC hints that DOC concentration and diffusion influences respiration in frozen soils.Tucker (2014) explicitly modeled DOC diffusion in the thin water films and consequently better reproduced observed respiration at temperatures just below freezing.Improving the simulated respiration at temperatures between −5 and 0 • C requires both the representation of the thawed organic matter and the effects of DOC diffusion.Incorporating a DOC carbon pool and DOC diffusion into SiBCASA is not within the scope of our current study and left as a future improvement to the frozen soil biogeochemistry parameterization.
For temperatures below −5 • C, the simulated respiration for the organic soil agrees with the incubation data, while the simulated respiration for the mineral soil is higher than observed.Our frozen biogeochemistry parameterization may not include some key processes that influence respiration below −5 • C. Eberling and Brandt (2003) found that the frozen samples trapped 10 % of the CO 2 produced, but this would apply to all temperatures and would not explain the higher simulated values below −5 • C. Intracellular desiccation or similar internal processes may slow down microbial activity and reduce respiration (Mikan et al., 2002).However, nearly all frozen carbon lies in the top 3 m of permafrost (Tarnocai et al., 2009;Hugelius et al., 2013).Since permafrost soils at these depths spend only a portion of the year at temperatures below −5 • C, the overall effect of this bias is small.
Whether a model needs the frozen biogeochemistry parameterization to represent substrate availability or whether the original Q 10f formulation would suffice depends upon the model application.If the model application focuses on short-term carbon fluxes on a seasonal to decadal timescale, the original Q 10f formulation would suffice.In such shortterm applications, the Q 10f formulation will produce realistic respiration fluxes, especially in winter.However, the Q 10f formulation does not work well in model applications focusing on 500-to 10 000-year timescales to study the buildup of frozen carbon or potential future emissions from thawing permafrost.For such long timescales, the model will need to account for substrate availability in order to correctly simulate the frozen permafrost carbon pools and correctly estimate future carbon fluxes from thawing permafrost.
When accounting for substrate availability, both the bulk and film frozen pools are required to build up or maintain the frozen permafrost carbon.A single set of frozen pools effectively reproduces the temperature effects represented by the original Q 10f formulation, but does not maintain the permafrost carbon pool.A single set of frozen pools suffers from numerical diffusion, which is an artificial dispersion of carbon in the model resulting from insufficient spatial finite differencing.For a single set of frozen pools, even though the permafrost was always below freezing, the simulated permafrost temperature and thus ϕ varied throughout the year such that carbon was artificially pumped from the frozen pools into the thawed pools.The amount of frozen carbon in the permafrost steadily decreased with effective turnover times of 200 to 500 years, the same result achieved using the original Q 10f formulation.Two sets of frozen carbon pools is the minimum number required to represent the physical separation of soil carbon located in the thin water films and minimizes artificial carbon loss due to numerical diffusion.The film frozen pools still have numerical diffusion, which can be decreased further by increasing the number of film frozen pools, but we leave such exploration as potential future work.Separating frozen carbon into film and bulk frozen pools is sufficient to minimize numerical diffusion and limit microbial activity to the substrate physically within thin water films.
Improving the frozen biogeochemistry parameterization will require additional measurements focusing on the effects of TOC and soil texture on respiration and VWC.The incubation studies we show here emphasized the lower temperature limits of microbial activity in frozen soil, with incubations at temperatures as low as −39 • C. Since the top 3 m of permafrost containing the bulk of the frozen carbon rarely reach such low temperatures, a more useful range for modelers would be 0 to −5 • C, where VWC and respiration show the greatest changes with temperature.Studies that quantify nutrient availability and diffusion in the thin water films would prove especially useful to quantify the non-linear response of respiration to VWC.For both incubation studies and VWC measurements, including ancillary data is very important in applying the data to models, especially the TOC, soil texture, soil bulk density, and water content.Also, the soil texture and organic content are much more important than where the sample was collected, so we need incubations and measurements that explore how TOC and soil texture influence VWC and respiration in frozen soils.Lastly, modelers need uncertainty estimates for the incubation and VWC measurements.The ideal performance for any model is to match observations within uncertainty, which indicate that model output and observations are statistically identical.This makes uncertainty estimates as important as the observations themselves, but none of the incubation and VWC studies shown here included uncertainty.

Conclusions
The frozen biogeochemistry parameterization successfully links soil physics, microbiology, and biogeochemistry to model substrate availability and associated effects on simulated respiration fluxes in frozen soils.The resulting simulated VWC is consistent with observed values and is strongly influenced by soil organic content.The simulated respiration fluxes as a function of temperature are consistent with observed values from incubation experiments and also depend very strongly on soil organic content.Future versions of the frozen biogeochemistry parameterization should account for additional, non-linear effects of substrate diffusion in thin water films on simulated respiration.Controlling respiration in frozen soils based on substrate availability allows us to maintain a realistic permafrost carbon pool by eliminating the continuous carbon loss seen in the original kinetic Q 10f formulation.The frozen biogeochemistry parameterization is a useful way to represent the effects of substrate availability on soil respiration in model applications that focus on century to millennial timescales in permafrost regions.
Figure1.A schematic of the new frozen biogeochemistry parameterization around a single, idealized soil particle.Soil carbon is divided into thawed carbon in the thin water films surrounding soil particles, film frozen carbon in the maximum extent of the thin liquid water film, and bulk frozen carbon in the pore spaces between soil particles.The thawed carbon represents available substrate for metabolism by microbes in the thin water film.
and b i is an empirical constant that depends on soil texture.Kurylyk and Watanabe (2013) compare several mathematical formulations for

Figure 2 .
Figure 2. The parameterized liquid water fraction, ϕ i , as a function of temperature for pure sand, silt, clay, and organic matter.The thin dashed lines show the corresponding values of critical water fraction, ϕ crit , defining the boundary between bulk freezing and thin film effects for each soil type.

Figure 4 .
Figure 4. Simulated daily average VWC for organic soil at 19 cm depth (a) and mineral soil at 54 cm depth (b) as a function of daily average temperature at Bonanza Creek.

Figure 5 .
Figure 5. Observed respiration as a function of temperature from the incubation studies in Table3.The Koppara and Plotnikovo incubation results appear as the same color.

Figure 6 .
Figure 6.Observed respiration from the incubation studies in Table 3 and simulated hourly respiration as a function of daily average temperature at Toolik Lake, Alaska, for organic soil (a) and mineral soil (b).

Figure 7 .
Figure 7. Simulated permafrost carbon in Northern Hemisphere permafrost with and without the frozen biogeochemistry parameterization.

Table 2 .
Transfer rules between thawed, frozen film, and frozen bulk pools.

Table 3 .
Summary of incubation data used to evaluate the frozen biogeochemistry parameterization.
* Numbers in bold were calculated from information supplied in the paper.theNorth Slope of Alaska (68.65 • N, 149.65 • W) and Bonanza Creek near Fairbanks, Alaska (64.80 • N, 148.00 • W).