Model responses to CO2 and warming are underestimated without explicit representation of Arctic small-mammal grazing.

We use a simple model of coupled carbon and nitrogen cycles in terrestrial ecosystems to examine how "explicitly representing grazers" vs. "having grazer effects implicitly aggregated in with other biogeochemical processes in the model" alters predicted responses to elevated carbon dioxide and warming. The aggregated approach can affect model predictions because grazer-mediated processes can respond differently to changes in climate compared with the processes with which they are typically aggregated. We use small-mammal grazers in a tundra as an example and find that the typical three-to-four-year cycling frequency is too fast for the effects of cycle peaks and troughs to be fully manifested in the ecosystem biogeochemistry. We conclude that implicitly aggregating the effects of small-mammal grazers with other processes results in an underestimation of ecosystem response to climate change, relative to estimations in which the grazer effects are explicitly represented. The magnitude of this underestimation increases with grazer density. We therefore recommend that grazing effects be incorporated explicitly when applying models of ecosystem response to global change.

Introduction

Despite evidence that animals can influence ecosystem carbon (C) and nutrient cycles (Schmitz 2014), the explicit incorporation of animals into terrestrial biogeochemical models is rare (Metclafe and Olofsson 2015). To maintain mass balance in these models without explicit representation of animals, the effects of animals have to be implicitly aggregated into other biochemical processes through model calibration (e.g., animal respiration included with other heterotrophic respiration). However, animal‐mediated processes can behave differently from the processes with which they are aggregated. For example, combining microbial and mammal respiration into a single value for heterotrophic respiration can cause problems, because warming generally increases respiration in microbes, but can slow respiration in mammals if the warming reduces the energy needed to maintain body temperature (Batzli et al. 1980). Here we examine the effects of aggregating grazer‐mediated processes in with other biogeochemical processes when modeling ecosystem response to elevated carbon dioxide (CO2) and warming. We use small‐mammal grazers in Arctic tundra as an example. However, the principles are general, and we conclude with a discussion of how our results might apply more generally to other grazers and other ecosystems.

Recent studies suggest that animals influence the response of tundra to climate change (Tuomi et al. 2019, Petit Bon et al. 2020). Experimental manipulations conducted across a range of tundra ecosystems have shown that, while warming or fertilization typically enhances above ground productivity and nutrient cycling in tundra, the presence of herbivores – including rodents, geese, and ungulates – can dampen or negate this response (e.g., Grellman 2002, Post and Pedersen 2008, Sjögersten and van der Wal 2008, Rinnan and Stark 2009, Cahoon et al. 2011, Kaarlejärvi and Hoset 2015, Leffler et al. 2019) or might enhance productivity (Gough et al. 2012). Furthermore, observational studies have shown that trophic interactions on the tundra strengthen under warmer conditions (McKinnon et al. 2010, Legagneux 2012), suggesting that animal influences on C cycling might be stronger in the future.

In other terrestrial ecosystems, animals are known to affect C and nutrient cycling (McNaughton 1985, McLaren and Jefferies 2004, Wilkinson and Sherratt 2016). The direct effects of animals vary among ecosystems, type of herbivore, and plant growth form (Jai et al. 2018) but have historically been thought of as small, relative to plant and microbial processes (e.g., Hairston and Smith 1960). Nevertheless, animals can accelerate nutrient cycles and influence plants and microbes indirectly by mediating chemical and biological processes and altering community structure and can thereby have a large influence on ecosystem C and nutrient processing (Pastor et al. 1988, Wardle et al. 2004, Zimov et al. 2009, Metcalfe 2014, Schmitz et al. 2014).

Although herbivore–vegetation models have been made for other ecosystems (e.g., Seagle and McNaughton 1993, Bennett 2003), we are aware of only one vegetation‐dynamics model – ArcVeg – that explicitly addresses the effect of an Arctic herbivore (caribou) on tundra biogeochemistry (Yu et al. 2011). This model indicates that grazing dampens the increase in plant biomass expected from warming soils and the consequent increase in nutrient cycling (Yu et al. 2011). These results suggest that the explicit inclusion of grazers in biogeochemical models could be necessary for predicting tundra responses to climate change.

Other modeling studies have addressed Arctic biogeochemical responses to climate change, but without the explicit representation of the effects of animals as separate from other C and nutrient‐cycling processes. These biogeochemical models indicate significant long‐term impacts of elevated CO2 and warming, but the model predictions differ on how these responses will ultimately affect net C source vs. sink activity. This source–sink disparity is due to uncertainty in the balance between elevated autotrophic and heterotrophic respiration (source) resulting from warming vs. enhanced photosynthesis (sink) resulting from the direct effects of elevated CO2 and warming on production and from the acceleration of nutrient cycles by warming (McKane et al. 1997, Rastetter and Ågren 1997, McGuire 2012, Pearce et al. 2015, Jiang et al. 2017). This trade‐off between source vs. sink activity is likely to be confounded by Arctic herbivores.

From the perspective of ecosystem biogeochemistry, aggregating herbivore effects in with other processes can be justified, because grazers perform several nutrient‐cycling processes that parallel other plant and microbial processes within ecosystems. Here we examine the aggregation effects for four such processes in relation to the response of ecosystems to elevated CO2 and warming:

Grazing mediates the transfer of plant C to detritus and soil organic matter (soil), and thereby acts in parallel with tissue senescence and litter fall.

Similarly, grazing transfers plant N to soil organic matter in parallel with tissue senescence and litter fall, but does so before the plants can resorb N.

Consumption of plant material and subsequent heterotrophic respiration by grazers parallels litter fall and the subsequent heterotrophic respiration resulting from processing of soil organic matter by microbes and other detritivores.

Finally, metabolic processing of plant matter consumed by grazers produces dissolved labile N in urine in parallel with litter fall and microbially mediated mineralization.

Even though these processes act in parallel, the grazer‐mediated and non‐grazer‐mediated processes might respond differently to climate change, or even in opposite directions. Furthermore, the cyclic dynamics of small‐mammal grazers in the Arctic might complicate the relative contributions of these parallel processes to ecosystem responses to elevated CO2 and climate change. Based on the modeling analysis we present below:

We hypothesize that aggregating the effects of small‐mammal grazers with other C and N cycling processes results in an underestimation of tundra responses to elevated CO2 and warming. For our model, after 100 years the underestimation of C sequestration in tundra ecosystems in response to elevated CO2 and warming is 50–80% relative to estimations in which the grazer effects are explicitly represented.

We hypothesize that although three‐to‐four‐year cycles in the density of small‐mammal grazers have measurable short‐term effects of tundra biogeochemistry (e.g., Olofsson and Tømmervik 2012), densities averaged over the grazer cycles can be used to assess long‐term responses of tundra to elevated CO2 and warming.

We use a simple model of coupled C and N cycles in ecosystems applied to the effects of small‐mammal grazers on the responses of moist acidic tundra to elevated CO2 and warming. Most of the data we use are for lemmings and voles, but the model applies to generic small‐mammal grazers in the Arctic, which we refer to as “voles” for simplicity. We apply the model both with vole densities explicitly represented and with vole densities unspecified, but their effects implicitly subsumed into other biogeochemical processes. In all applications of the model, we assume voles are present on the landscape. The model applications differ only in the way that vole and other biogeochemical processes are separated from one another.

Methods

Model

We use a model developed by Rastetter et al. (2020) to examine recovery of ecosystems from disturbances that remove vegetation (Box 1; Rastetter et al. 2021). We have modified that model to account for temperature sensitivity of six metabolic processes (photosynthesis, autotrophic and heterotrophic respiration, plant and microbial N uptake, and N mineralization). We have also modified it to account for the effects of voles on the transfer of C and N from vegetation to soil organic matter and the transfer of N in urine from vegetation to inorganic N (although not all N in urine is inorganic, it is labile, and we treat it as inorganic). The basic model is fully described in Rastetter et al. (2020); here we describe only the changes to that model for the current analyses.

Variables and parameters defined in Table 1.

Table 1

Model variables and parameters.

	Symbol	Value	Units
C and N stocks
Vegetation C	B C	878	g C/m2
Detritus and soil organic C	D C	19,452	g C/m2
Vegetation N	B N	20.6	g N/m2
Detritus and soil organic N	D N	831	g N/m2
Inorganic N	N	0.27	g N/m2
Processes and constraints
Allometric constraint	S	243.75	g C/m2
Vegetation stoichiometric constraint	Ψ	1	None
Soil stoichiometric constraint	Φ	1	None
Photosynthesis	P s	430	g C·m−2·yr−1
Autotrophic respiration	R a	215	g C·m−2·yr−1
Litter‐fall C	L itC	215	g C·m−2·yr−1
Heterotrophic respiration (excluding voles)	R h	213.07	g C·m−2·yr−1
Vegetation N uptake	U N	5.3800	g N·m−2·yr−1
Litter‐fall N	L itN	5.3800	g N·m−2·yr−1
Gross N mineralization	N min	19.7310	g N·m−2·yr−1
N immobilization	U Nm	14.4824	g N·m−2·yr−1
Inorganic N losses	Q DIN	0.0016	g N·m−2·yr−1
Refractory N losses	Q NR	0.1314	g N·m−2·yr−1
Refractory C losses	Q CR	1.93	g C·m−2·yr−1
C removed from vegetation by voles	G C	0	g C·m−2·yr−1
Vole respiration	R V	0	g C·m−2·yr−1
C added to soil organic matter by voles	L VC	0	g C·m−2·yr−1
N removed from vegetation by voles	G N	0	g N·m−2·yr−1
C added to soil organic matter by voles	L VN	0	g N·m−2·yr−1
Vole N transfer vegetation to inorganic soil	V UN	0	g N·m−2·yr−1
Driver variables
Atmospheric CO2	C a	400	μmol/mol
Temperature	T	10	°C
N inputs	N in	0.1330	g N·m−2·yr−1
Voles	V	0	voles/ha
Parameters
Allometric parameter 1	α	0.002231	m2·g−1 C
Allometric parameter 2	γ	0.01100	m2·g−1 C
Optimum vegetation C:N	q B	42.62	g C·g−1 N
Optimum soil C:N	q D	23.41	g C·g−1 N
Photosynthesis rate parameter	g C	1.423	yr−1
CO2 half‐saturation constant	k C	100.0	μmol/mol
Photosynthesis Q‐10	Q 10Ps	1.550	None
Autotrophic respiration constant	r B	0.09069	yr−1
Autotrophic respiration Q‐10	Q 10Ra	2.700	None
Vegetation C turnover rate constant	m CB	0.2449	yr−1
Vegetation N‐uptake rate parameter	g N	0.05191	g N·g−1 C·yr−1
Vegetation N half‐saturation constant	k N	1.000	g N/m2
Vegetation N‐uptake Q‐10	Q 10U	2.000	None
Vegetation N turnover rate constant	m NB	0.2612	yr−1
Heterotrophic respiration constant	r D	0.003651	yr−1
Heterotrophic respiration Q‐10	Q 10Rh	3.000	None
Microbial N‐uptake rate parameter	g Nm	0.001796	g N·g−1 C·yr−1
Microbial N half‐saturation constant	k Nm	1.000	g N/m2
Microbial N‐uptake Q‐10	Q 10Um	1.950	None
Vole nesting material	n V	22	g C·vole−1·yr−1
Vole C ingestion rate	g V	3,512	g C·vole−1·yr−1
Temperature slope vole metabolism	εV	52	g C·vole−1·°C−1·yr−1
Summer to annual temperature correction	T 0	10	°C
C:N of vole forage and nest material	q V	19.15	g C·g−1 N
Vole base respiration rate	r V	0.3	None
Vole urine N production rate	m NV	11.00	g N·vole−1·yr–1
Soil organic N turnover constant	m Nm	0.01099	yr–1
Soil organic N turnover Q‐10	Q 10m	2.160	None
C:N of DOM loss	q DOM	14.69	g C·g−1 N
N loss‐rate parameter	βN	0.005926	yr−1
Refractory N loss parameter	βNR	0.0001581	yr−1

Variable values are for the initial steady state with the aggregated representation of vole effects. Ψ and Φ are assumed to equal 1 under this steady state. Q 10 values are as reported in the main text. Other values are from Pearce et al. (2015) or are fitted to analogous functions in Pearce et al. (2015). Parameters are listed to four significant digits. DOM, dissolved organic matter.

Temperature response of metabolic processes

Because we use an annual time step in the model (i.e., no seasonality) and restrict warming to 5°C above current temperatures, we use a simple Q 10 function to simulate temperature responses rather than more complex formulations (e.g., Heskel 2016 or Carey 2016). We have therefore modified the photosynthesis, autotrophic respiration, heterotrophic respiration, plant N uptake, microbial N uptake, and N mineralization in the Rastetter et al. (2020) model to increase exponentially with warming (Box 1: Eqs. 9, 10, 13, 14, 21, 22).

Vole grazing

We drive the model by specifying vole density in each year (V). Consistent with values and cycle frequency reported in the literature (Batzli et al. 1980, Krebs and Boonstra 1995, Korpimaki et al. 2004, Pitelka and Batzli 2007, Krebs 2013, Ehrich 2020), we use a randomly generated time series of vole abundance with peaks every three or four years, with abundances at the peak ranging from 90 to 110 voles/ha, minimum abundances ranging from 8 to 12 voles/ha, and a mean vole abundance over the full time series of 40 voles/ha (Fig. 1). We chose to take this prescribed approach to vole density because the drivers of vole cycles are not fully understood (Korpimaki et al. 2004, Prevedello et al. 2013, Oli 2019) and are likely to include a top‐down component (Pitelka and Tomich 1955, Hairston et al. 1960, Krebs 2013), which is well beyond the domain of our model. For convenience, we specify vole density in voles/ha and correct to m−2 units by dividing by 10,000 m2/ha (Eqs. 15, 17, 18).

Fig. 1

Vole cycle used to simulate ecosystem response to grazing. Vole abundance is randomly generated with peaks every three or four years, with abundances at the peak ranging from 90 to 110 voles/ha, minimum abundances ranging from 8 to 12 voles/ha, and a mean vole abundance of 40 voles/ha. Upper panel shows the first 30 years of the time series. Bottom panel is the full 200‐year time series.

The removal of C from vegetation for nest building and ingestion are lumped into a single process for our model (Eq. 15). We assume this C removal is proportional to the specified vole density but decrease the per capita rate of ingestion with warming to account for decreased energy requirements to maintain body temperature (Eq. 15; Batzli et al. 1980). Vole respiration is proportional to the ingestion component of C removal from vegetation by voles and therefore also declines with warming (Eq. 17). We do not account for other temperature responses such as those associated with cold or heat stress. We assume a constant C:N ratio of material removed from vegetation by grazers (Eq. 16) but, because of the respiration loss of C and the N transferred to inorganic N in urine, the C:N of material removed from vegetation and the C:N of material added to the soil organic matter differ. Finally, we assume urine losses of N are proportional to vole density (Eq. 18).

Model calibration

We calibrate the model to be consistent with the C and N stocks and process rates compiled by Pearce et al. (2015) for the Multiple Element Limitation (MEL) model (Table 1). Because voles are part of the ecosystem, the effects of voles on tundra C and N stocks and fluxes are implicitly included in the data compiled by Pearce et al. (2015). Therefore, by using these data to calibrate (fit) the model without explicit vole representation, we are implicitly aggregating those vole effects in with the parallel ecosystem processes described above. In the calibrations in which voles are explicitly represented, we assume a constant vole density and that the ecosystem is in steady state. We then specify the vole effects directly and subtract these effects from the parallel ecosystem processes before calibration (Table 2). The combined rates of vole‐mediated processes plus the parallel ecosystem processes are therefore identical in all calibrations (rows labeled “Total” in Table 2), therefore providing the basis of comparison for assessing the consequences of explicit vs. aggregated representation of vole grazing.

Table 2

Variable and parameter changes to accommodate the effect of explicit representation of voles.

Symbol	Total or PAR	Vole effects aggregated in with other processes	Explicit vole representation with 40 voles/ha	Explicit vole representation with 100 voles/ha	Units
L itC		215	200.864 (−6.6%)	179.66 (−16.4%)	g C·m−2·yr−1
G C		0	14.136	35.34	g C·m−2·yr−1
	Total	215	215	215	g C·m−2·yr−1
m CB	PAR	0.2449	0.2288 (−6.6%)	0.2046 (−16.4%)	yr−1
R h		213.07	208.8556 (−2.0%)	202.534 (−4.9%)	g C·m−2·yr−1
R v		0	4.2144	10.536	g C·m−2·yr−1
	Total	213.07	213.07	213.07	g C·m−2·yr−1
r D	PAR	0.003651	0.003658 (−2.0%)	0.003471 (−4.9%)	yr−1
L itN		5.38	4.6418 (−13.7%)	3.5346 (−34.3%)	g N·m−2·yr−1
G N		0	0.7382	1.8454	g N·m−2·yr−1
	Total	5.38	5.38	5.38	g N·m−2·yr−1
m NB	PAR	0.2612	0.2253 (−13.7%)	0.1716 (−34.3%)	yr−1
V UN		0	0.044	0.110	g N·m−2·yr−1
N min		19.731	19.687 (−0.2%)	19.621 (−0.6%)	g N·m−2·yr−1
	Total	19.731	19.731	19.731	g N·m−2·yr−1
m Nm	PAR	0.01099	0.01097 (−0.2%)	0.01093 (−0.6%)	yr−1

Values in parentheses are the percent change from the values used in the implicit‐vole representation with vole effects aggregated in with parallel ecosystem processes. “Total” is the total of the vole‐mediated and the parallel ecosystem process in the two preceding rows. “PAR” is the parameter in the equation for the parallel process in the preceding rows that was modified to accommodate explicit representation of voles.

In the calibrations, first we set allometric, C:N, Q 10, half‐saturation, and vole‐related parameters (derivation of these parameter values is presented in Appendix S1). We then set the rate parameters for each process so that flux rates are consistent with rates reported in Pearce et al. (2015). Because annual rates of plant and microbial processes are dominated by growing‐season rates, we use average summer temperature (10°C) to calibrate the model; in any case, because of the Q 10 formulation, once calibrated to a specified temperature, model responses are sensitive to changes in temperature, not to the temperature value itself. For vole processes, we correct this summer temperature to average annual temperature with an off‐set (Eq. 15).

We made three calibrations (Tables 3, 4). In calibration I, vole densities are not explicitly specified, and we assume that vole‐mediated processes can be implicitly represented by aggregating them with the parallel biogeochemical processes described in the introduction above (Table 2). In this calibration, we therefore set the number of voles (V) in the model to zero but incorporate the effects of voles in with the parallel processes through the calibration. In calibration II, we set the number of voles to 40 voles/ha so that vole‐mediated processes are explicitly represented, and the number of voles is the mean abundance of voles we use in our simulated vole cycle (described above). In calibration III, we set the number of voles to 100 voles/ha, which is the mean peak vole abundance in our simulated vole cycle. In calibration III, the parallel ecosystem process rates are decreased proportionally more than in calibration II to account for the higher vole density (Table 2).

Table 3

Calibrations.

Calibration	Vole representation	Vole density
I	Vole effects aggregated in with other processes	Unspecified vole density, but vole effects subsumed into litter‐fall C and N, heterotrophic respiration, and N mineralization in the calibration (Table 2)
II	Explicit vole representation	40 voles/ha
III	Explicit vole representation	100 voles/ha

Table 4

Simulations.

Simulation	Calibration	Description	Figure
Set 1
1	II	Constant 40 voles/ha	Fig. 2 dotted black lines
2	II	Voles cycling as in Fig. 1	Fig. 2 cycling black solid lines
3	II	Voles cycling as in Fig. 1 for 10 yr then held constant at 100 voles/ha	Fig. 2 dashed red lines
4	II	Voles cycling as in Fig. 1 for 10 yr then held constant at 0 voles/ha	Fig. 2 dash‐dotted blue lines
Set 2
5	I	Vole density unspecified, linear increase of CO2 from 400 to 800 μmol/mol over 100 yr	Figs. 4, 5, 6 dotted line, left column
6	I	Vole density unspecified, linear increase in temperature from 10 to 15°C over 100 yr	Figs. 4, 5, 6 dotted line, middle column
7	I	Vole density unspecified, linear increase of CO2 from 400 to 800 μmol/mol and temperature from 10 to 15°C over 100 yr	Figs. 4, 5, 6 dotted line, right column
8	II	Vole density cycling as in Fig. 1, linear increase of CO2 from 400 to 800 μmol/mol over 100 yr	Figs. 4, 5, 6 blue solid line, left column
9	II	Vole density cycling as in Fig. 1, linear increase in temperature from 10 to 15°C over 100 yr	Figs. 4, 5, 6 blue solid line, middle column
10	II	Vole density cycling as in Fig. 1, linear increase of CO2 from 400 to 800 μmol/mol and temperature from 10 to 15°C over 100 yr	Figs. 4, 5, 6 blue solid line, right column
11	III	Constant 100 voles/ha, linear increase of CO2 from 400 to 800 μmol/mol over 100 yr	Figs. 4, 5, 6 dashed red line, left column
12	III	Constant 100 voles/ha, linear increase in temperature from 10 to 15°C over 100 yr	Figs. 4, 5, 6 dashed red line, middle column
13	III	Constant 100 voles/ha, linear increase of CO2 from 400 to 800 μmol/mol and temperature from 10 to 15°C over 100 yr	Figs. 4, 5, 6 dashed red line, right column

All but four of the parameters have the same values in all three calibrations. To maintain the same steady state in all three calibrations, we adjust the values of these four parameters to compensate for how voles are represented in the model (m CB, m NB, r D, and m Nm; rows labeled “PAR” in Table 2). These parameters are adjusted so that the rates of C and N litter losses, heterotrophic soil respiration, and gross N mineralization all decrease to compensate for the parallel vole‐mediated C and N fluxes in calibrations II and III in which voles are explicitly represented. Because we calibrate to the same data set (Pearce et al. 2015), the overall C and N stocks and cycling rates are identical for these three calibrations (rows labeled “Total” in Table 2).

Simulations

We run a total of 13 simulations in two sets (Tables 3, 4; Rastetter et al. 2021, ).

In the first set of simulations, we assume that the average vole density is 40 voles/ha and use calibration II with vole effects explicitly represented (Table 3). We then run four 200‐year simulations with no change in either CO2 or temperature. We drive the model with: (1) voles held constant at 40 voles/ha; (2) voles cycling on the three‐to‐four‐year cycle between 8 and 110 voles/ha for 200 years; (3) voles cycling for 10 years, followed by maintenance of a constant vole density of 100 voles/ha (equivalent to adding voles to the ecosystem); and (4) voles cycling for 10 years, followed by maintenance of a constant vole density of 0 voles/ha (equivalent to removing voles from the ecosystem).

This first set of simulations serves two purposes. First, it illustrates the effects of long‐term changes in vole density and thereby draws the distinction between adding or removing voles from calibrating the model assuming high or low vole density. Second, it allows us to assess the potential long‐term effects of voles if their numbers were maintained at high or low levels. We can thereby address the question: “Do the simulated changes in the ecosystem approach their potential changes during peaks and troughs in the vole cycle?”

The second set of simulations is to address our central question about aggregated vs. explicit representations of grazer effects on ecosystem responses to climate change (Table 4). We run nine 100‐year simulations in a two‐factor design. The first factor relates to how vole effects are represented in the three calibrations (Table 3) and vole abundance:

Calibration I (aggregated) and vole abundance subsumed in the calibration of the parallel processes and therefore assumed constant but unspecified (although V is set to 0 in the model, vole effects are aggregated in with the parallel ecosystem processes).

Calibration II (40 voles/ha) and vole abundance cycling on the three‐to‐four‐year cycle between 8 and 110 voles/ha.

Calibration III (100 voles/ha) and vole abundance held constant at 100 voles/ha.

The second factor relates to climate change:

A linear increase in atmospheric CO2 from 400 to 800 μmol/mol over 100 years.

A linear increase in temperature from 10 to 15°C over 100 years.

A linear increase in both atmospheric CO2 from 400 to 800 μmol/mol and temperature from 10 to 15°C over 100 years.

Results

Set 1: Effects of vole cycling and adding or removing voles

Simulation 1: Effects of holding voles constant at the calibration abundance

Because the model was calibrated to a steady state with 40 voles/ha, all ecosystem C and N stocks and fluxes remained constant when vole abundance was held at 40 voles/ha in the 200‐year simulations (dotted horizontal lines in Fig. 2). This simulation only serves to illustrate the stability of the model and to serve as a control to which the other simulations can be compared.

Fig. 2

Simulated changes in plant and soil organic carbon (C) and nitrogen (N) in response to constant vole abundance, vole cycling, vole density maintained at 100 voles/ha, and vole removal (see Table 4). The thin dotted black lines are the steady‐state values if vole density is held at 40 voles/ha (simulation 1). Solid black lines are the responses to the vole cycle depicted in Fig. 1 (simulation 2). Dashed red lines are the responses to the same vole cycle and then vole density maintained at 100 voles/ha after year 10 (simulation 3). Dashed‐dotted blue lines are the responses to the same vole cycle and then removal of voles after year 10 (simulation 4).

Simulation 2: Effects of vole cycling on plant and soil C and N

In the 200‐year simulations with vole abundance cycling, the plant and soil C and N stocks cycle at the same three‐to‐four‐year frequency as the voles (Fig. 2). In addition, there are some longer term dynamics in these stocks associated with the autocorrelated nature of plant production and the legacy of the random variations in the vole cycle. Despite these dynamics, vole cycling does not cause the plant and soil C and N stocks to diverge far from the values to which they are calibrated (dotted and solid lines in Fig. 2).

The plant biomass cycles out of phase with the vole cycle. The lowest plant C and N values occur in years of peak vole numbers and the highest plant C and N values occur three or four years after peak vole numbers or the year prior to the subsequent vole peak (Fig. 3). This phase shift in the plant relative to vole cycles as well as the magnitude of the plant C cycle (20–30 g C/m2 peak to trough) are roughly consistent with the phase and magnitude of the cycles reported by Olofsson et al. (2012). In addition, the dependence of plant production on plant biomass results in a strong autocorrelation in the plant C and N time series, which in turn results in longer term dynamics less clearly tied to the vole cycle (Fig. 2). This autocorrelation is reflected in the strong positive correlation between the plant C and N and their respective values at the time of the previous peak in vole numbers (open dots remain high and closed dots remain low in Fig. 3).

Fig. 3

Plant carbon (C) and nitrogen (N) recovery following peak vole abundance in simulation 2. The plant C and N of 875 g C/m2 and 20.2 g N/m2 were selected to partition the recovery time series into two approximately equal‐sized groups based on their values at the time of the previous vole peak (time 0 on x axis). The levels of C and N during this recovery depend not only on peak vole abundance, but also on the degree of recovery following the previous vole cycle. Because the biomass consumed is proportional to vole abundance and not to plant biomass, if plants recover to a higher level following the previous cycle (white dots), then they begin and maintain recovery in the current cycle at a higher level relative to plants that recovered to a lower level in the previous cycle (black dots). This autocorrelation results in the longer term dynamics in Fig. 2 for simulation 2 in which vole abundance cycled. The recovery in any cycle also depends on vole abundance and the duration of the vole cycle (higher plant recovery in a four‐year cycle than a three‐year cycle).

The dynamics of soil C and N stocks are closely tied to the plant dynamics. Because N inputs to the ecosystem are less than 3% of the annual plant requirement (Table 1), plant recovery from vole outbreaks relies almost exclusively on N from soil organic matter. As a consequence, the three‐to‐four‐year soil N cycles are directly out of phase with plant N cycles and the longer term dynamics are also opposite those in plant N (Fig. 2). Soil C also cycles out of phase with plant C, but the relation is not as strong as it is for N. However, because soil C is ultimately derived from plant C, the longer term dynamics in plant C are paralleled in the soil C following approximately a 9‐year lag (Fig. 2).

Simulations 3 and 4: Effects of removing or adding voles

When voles are removed from the ecosystem, plant C and N increase by approximately 13%, or increase by 116 g C/m2 and 2.7 g N/m2. Because of the reliance of plants on soil N, soil N deceases by almost the same absolute amount as the plants gain, 2.5 g N/m2. However, the amount of N in the soil is so large that this loss amounts to only an approximately 0.3% loss. The increase in plant biomass results in higher litter inputs to the soil. Consequently, soil C increases by approximately 4% or 727 g C/m2. The gain of soil C and loss of soil N widens the soil C:N ratio by approximately 4%, which in turn increases microbial N immobilization into soil organic matter (effect of Φ in Eq. 14).

When voles are increased and then held constant at 100 voles/ha, plant C and N decrease by approximately 20%, or 175 g C/m2 and 4 g N/m2. Again, because of the tight cycling of N in the ecosystem, soil N increases by almost the same absolute amount as the plants lose, 3.7 g N/m2 (0.4%). The loss of plant biomass translates into lower litter inputs to soil and a large absolute decrease in soil C, 1,170 g C/m2 (6%). Because of the increase in soil N and decrease in soil C, the soil C:N narrows by approximately 6%, which in turn decreases microbial N immobilization into soil organic matter (effect of Φ in Eq. 14).

In the simulations in which vole density is increased or decreased and then held constant, it takes 10–20 years for the vole effects to reach their largest deviation from the steady state and another 60–90 years for those effects to stabilize. Because of this long response time, the potential effects of voles on tundra biogeochemistry cannot be approached if vole abundance cycles on a three‐to‐four‐year cycle. Indeed, when voles are cycling, the magnitude of these effects relative to the peak effects of the long‐term increase or decrease in vole abundance is only approximately 12% for plant C, 2% for soil C, and 20% for both plant and soil N.

Set 2: Effects of aggregated vs. explicit representations of vole effects

Simulations 5, 8, and 11: Responses to increasing CO2

Predicted responses to increased CO2 do not differ substantially between the aggregated model in which vole effects are implicitly aggregated in with other biogeochemical processes through model calibration (calibration I) and the distributed model in which vole effects are explicitly represented (calibrations II and III: Figs. 4, 5, 6). The only process affected by elevated CO2 is photosynthesis. However, because the plants are strongly N limited, elevated CO2 increases net primary production (NPP) by only 10–11% in both aggregated and distributed simulations (Fig. 6). This increase in production translates into an approximately 11–12% increase in biomass, again in both aggregated and distributed simulations (Fig. 4). The increase in production results in only about a 4% increase in soil C in all the simulations. This increase in soil C is a small relative change but, because soil has such a large fraction of the organic matter, it is a large absolute change amounting to about 90% of the total change in ecosystem C.

Fig. 4

Simulated changes in plant, soil, and total ecosystem C with a linear increase in CO2 from 400 to 800 μmol/mol over 100 years, a linear warming from 10 to 15°C over 100 years, and both a linear increase in CO2 from 400 to 800 μmol/mol and a linear warming from 10 to 15°C over 100 years (see Table 4). Different trajectories indicate responses with vole effects aggregated with other biogeochemical processes (dotted black lines; simulations 5, 6, and 7), voles cycling between 8 and 110 voles/ha on a three‐to‐four‐year cycle (solid blue lines; simulations 8, 9, and 10), and a constant 100 voles/ha (dashed red lines; simulations 11, 12, and 13).

Fig. 5

Simulated changes in plant, soil, and total ecosystem N with a linear increase in CO2 from 400 to 800 μmol/mol over 100 years, a linear warming from 10 to 15°C over 100 years, and both a linear increase in CO2 from 400 to 800 μmol/mol and a linear warming from 10 to 15°C over 100 years (see Table 4). Different trajectories indicate responses with vole effects aggregated with other biogeochemical processes (dotted black lines; simulations 5, 6, and 7), voles cycling between 8 and 110 voles/ha on a three‐to‐four‐year cycle (solid blue lines; simulations 8, 9, and 10), and a constant 100 voles/ha (dashed red lines; simulations 11, 12, and 13).

Fig. 6

Simulated changes in gross primary (GPP), net primary (NPP), and net ecosystem production (NEP) with a linear increase in CO2 from 400 to 800 μmol/mol over 100 years, a linear warming from 10 to 15°C over 100 years, and both a linear increase in CO2 from 400 to 800 μmol/mol and a linear warming from 10 to 15°C over 100 years (see Table 4). Different trajectories indicate responses with vole effects aggregated with other biogeochemical processes (dotted black lines; simulations 5, 6, and 7), voles cycling between 8 and 110 voles/ha on a three‐to‐four‐year cycle (solid blue lines; simulations 8, 9, and 10), and a constant 100 voles/ha (dashed red lines; simulations 11, 12, and 13).

The amount of N entering the ecosystem is too small to support even the small gain in plant C in response to elevated CO2. The gain is instead supported by a net transfer of 0.9–1.3 g N/m2 from soil to plants over the 100‐year simulations. The amount of N transferred from soil to plants is about the same in all three simulations (Fig. 5). Elevated CO2 increases the C:N ratio of the plants, which in turn increases N uptake (effect of Ψ in Eq. 10). However, the increase in soil C:N resulting from increased litter inputs also increases microbial N uptake (effect of Φ in Eq. 14). This competition between plants and microbes for N limits the ecosystem response to elevated CO2. Again, the effects of aggregated vs. explicit representation of voles on this response to elevated CO2 are negligible (Figs. 4, 5, 6).

Simulations 6, 9, and 12: Responses to warming

In contrast, the effects of aggregated vs. explicit representation of voles on the response to warming are large (Figs. 4, 5, 6). Warming not only stimulates photosynthesis, it also stimulates autotrophic and heterotrophic respiration, and, more importantly, it stimulates the N cycle in three places: (1) N mineralization, (2) microbial N uptake, and (3) plant N uptake. A major effect of this stimulation of the N cycle is an increase in net N mineralization, a resulting relaxation of N limitation on plant growth, and a large increase in plant biomass. The increase in plant production increases litter inputs to soils, which in turn mitigates soil C losses. In addition, the increased production allows the ecosystem to accumulate a small amount of N (<0.3 g N/m2; Fig. 6). The effects of this chain of events are much stronger in the simulations in which voles are explicitly represented than in the simulations with the aggregated model and the effects are stronger when the model is calibrated assuming higher vole densities (response stronger for calibration III [100 vole/ha] than for calibration II [40 voles/ha]). Therefore, explicit representation of vole effects results in an amplification of the predicted transfer of N from soil to plants, larger predicted gains in plant C or higher predicted retention of soil C, and higher predicted rates of gross primary production (GPP), NPP, and net ecosystem production (NEP). If the simulations are run with vole density held constant at 40 voles/ha, the C and N stocks and fluxes follow the same general patterns as in the simulations with the three‐to‐four‐year vole cycle (data not shown). The size of the differences between the simulations with constant 40 voles/ha and cycling vole density are about the same as those of the cycle simulation from the steady state with no climate change (Fig. 2). In addition, the temperature effects on vole consumption and respiration (Box 1, Eqs. 15, 17) have only a small effect on this general pattern (simulations rerun with εV = 0 resulted in <2% difference in C and N stocks; data not shown).

Simulations 7, 10, and 13: Responses to increasing CO2 and warming

The effects of elevated CO2 and warming are slightly amplified when the two are combined (the two interact synergistically). Under both elevated CO2 and warming, GPP is 12% and NPP is 8% higher than the sum of the changes in GPP and NPP under elevated CO2 alone and warming alone (Fig. 6). The net transfer of N from soil to plants is about 3% stronger and the increase in plant C is about 7.5% stronger (Figs. 4, 5). Overall, because the response to CO2 alone is so much smaller that the response to warming alone, the response to the two combined is dominated by the warming response. The consequences of aggregated vs. explicit representations of vole effects are therefore the same as in the warming simulations.

In our analysis we assume vole density is top‐down controlled and therefore does not increase with plant production. However, if the average vole density during the cycle increases in proportion to NPP (~80% over 100 year), some of the increased production with elevated CO2 and warming is consumed and the increase in plant biomass is about 7.6% lower than when the average vole density remains constant (data not shown). The increase in vole density causes soil C to decrease by 0.9% rather than increase by 0.4%.

Discussion

Our analysis indicates that failure to explicitly represent small‐mammal grazers (voles) in biogeochemical models can result in an underestimation of the response of Arctic ecosystems to climate warming but has only a small effect on the response to elevated CO2 (Figs. 4, 5, 6). Underestimation of the warming response increases with the assumed density of voles used to calibrate the model. Although cycling of vole density has short‐term effects on ecosystem stocks and fluxes, it neither amplifies nor dampens the underestimation in the long‐term response to warming.

Why is the explicit representation so important?

Before addressing this question, we again emphasize the distinction between explicit representation of voles and adding voles to the ecosystem. Adding voles to the ecosystem accelerates nutrient cycling by increasing the transfer of nutrients from plants to soil. Such an acceleration of nutrient cycling might be expected to increase the responsiveness to elevated CO2 and warming. However, in our analysis of the response to elevated CO2 and warming we do not add voles, we simply change how the voles are represented in the model. Vole effects are either implicitly aggregated in with other processes or they are explicitly represented. In all three of our calibrations, the total amounts of C and N removed from vegetation, the total heterotrophic respiration, and the total mineralization of N are identical (Table 2). Therefore, the magnitudes of vole‐mediated processes plus the parallel ecosystem processes are represented identically in all three calibrations. Furthermore, in the analysis in which we did add voles, any acceleration of nutrient cycles by vole activity is transient; our analysis indicates that the net effect of adding voles is to transfer N from plants, with relatively high N turnover, to soil, with slower N turnover (Fig. 2); although it would be impossible to detect the 0.6% change in soil N predicted by our model. Therefore, the long‐term effect of adding voles is to slow the nutrient cycle. Furthermore, adding and maintaining 100 voles/ha resulted in a loss of over one kilogram of C from the ecosystem (Fig. 2). Therefore, the effect of adding grazers is to decrease ecosystem C whereas explicitly representing voles in the model is to increase the estimate of C gain with climate change (Fig. 4).

Why are the differences in responses to elevated CO2 so small between aggregated vs. explicit representation of vole effects?

Elevated CO2 stimulates only one ecosystem process, photosynthesis (Fig. 7). The associated increase in C gain increases biomass and leaf area, which further stimulates photosynthesis (Fig. 7: ↑C a → ↑P s → ↑B C → ↑S → ↑P s). However, there is a much stronger negative feedback associated with the change in stoichiometry; increased photosynthesis increases biomass C and consequently increases vegetation C:N, which feeds back to decrease photosynthesis (↑C a → ↑P s → ↑B C → ↑Ψ → ↓P s). Tissue and Oechell (1987) used this stoichiometric feedback to argue why tussock tundra exposed to elevated CO2 alone had only a transient increase in production. Without an increase in the N supply to vegetation, the CO2‐stimulation of photosynthesis cannot be maintained. However, the increase in vegetation C:N increases the litter‐fall C:N and consequently the soil C:N, which in turn decreases net N mineralization and the supply of N to plants (↑C a → ↑P s → ↑B C → ↑Ψ → ↓L itN → ↓D N → ↑Φ → ↓(N min − U Nm)). Because none of the steps in this chain were modified to incorporate voles explicitly in the model calibration, this N feedback is about the same for both explicit and aggregated representation of voles in the model and therefore does not have a large effect on the relative responses with and without explicit representation of vole effects.

Fig. 7

Causal‐chain diagram for the model in Box 1. Arrows indicate causal links: a red arrow marked with a “+” indicates that an increase in the variable at the tail of the arrow will cause an increase in the variable at the head of the arrow; a blue arrow marked with a “−” indicates that and increase in the variable at the tail of the arrow will cause a decrease in the variable at the head of the arrow. Symbols are defined in Table 1. Symbols in boxes are C and N stocks, symbols in circles are driver variables, and other symbols are either processes or allometric and stoichiometric constraints. The four causal links shown as dashed arrows are the links that were weakened in the calibration to accommodate vole‐mediated processes in the simulations with explicit representation of vole density (see Table 2). The temperature (T) and vole (V) drivers are shown three times to avoid overcomplicating the diagram.

Why are the predicted responses to warming stronger when vole effects are explicitly represented in the model than when they are aggregated with other processes?

In the model, warming stimulates six processes: photosynthesis, autotrophic and heterotrophic respiration, N mineralization, microbial N uptake, and plant N uptake (Fig. 7). Although warming decreased the energy cost of thermoregulation and therefore decreases vole consumption of plants (Eq. 15) and vole respiration (Eq. 17), we found that this effect is too small to explain the differences between simulations with vs. without voles explicitly represented (accounting for <2% of the response in C and N stocks).

Among these many effects of warming in the model, the main effect that results in the accumulation of plant C in simulations with both explicit and aggregated representations of vole effects is the release from N limitation through the stimulation of net N mineralization (Fig. 7: ↑T → ↑(N min − U Nm) → ↑N → ↑U N → ↑B N → ↓Ψ →↑P s). Therefore, one effect of this mobilization of soil N is for both C and N uptake by plants to increase, causing plant biomass to accumulate. Because N mineralization (N min) was decreased in calibrations II and III to accommodate the explicit representation of voles (Table 2), this warming‐induced growth in plant biomass is about 0.8% (40‐vole calibration II) to 2% (100‐vole calibration III) weaker in the simulations with the explicit representation of vole effects. These simulations nevertheless accumulate more, not less, biomass.

The main reason that plant C and N accumulation differed between simulations with aggregated vs. explicit representations of voles is the change made to litter‐fall rates to accommodate the voles in calibrations II and III (Table 2). Litter fall is not directly stimulated by warming (Eqs. 11, 12), but it does increase as plant biomass increases. This increase in litter fall in turn limits the amounts of C and N that can accumulate in plants. However, the fraction of plant C and N lost in litter fall was decreased in calibrations II and III in which voles are explicitly represented to accommodate the C and N consumed by voles (Table 2: m CB and m NB were decreased). Consumption by voles does not increase with plant biomass (Eq. 15), and vole density does not increase with plant biomass because we assume top‐down control on voles and therefore use a prescribed vole density. As a consequence, the increase in litter fall as biomass accumulates is smaller with voles explicitly represented than in the aggregate representation. The accumulation of C and N in vegetation is therefore larger with the explicit inclusion of vole effects than in the simulation in which vole effects are aggregated with other processes. When we do allow vole density to increase in proportion of the increase in NPP (~80% over 100 years), the increase in plant biomass is less than 8% lower because of consumption and the small increase in soil C (<0.5%) becomes a small decrease (<1%).

In addition, because the C:N ratio of forage (19.15 g C·g−1 N) is lower than the C:N ratio of litter (40 g C·g−1 N), the fraction of vegetation N lost in litter fall was decreased more than the fraction of vegetation C lost in litter fall in calibrations II and III with explicit vole representations (Table 2; e.g., 13.7% decrease in litter N vs. 6.6% decrease in litter C with 40 vole/ha). As a consequence, as vegetation biomass accumulates, the litter‐fall C:N ratio increases more with the explicit representation of voles than with aggregated representation of vole effects. Soil organic C therefore increases more with the explicit vole representation than without it (Fig. 4), and soil organic N decreases more with the explicit vole representation than without it (Fig. 5).

Why is there a synergistic response to elevated CO2 and warming in combination?

If the response to CO2 was stronger so that there was a substantial increase in plant biomass and plant C:N ratio, then the feedback associated with litter fall would have come into play and differences between explicit and aggregated representations of vole activity would have made more of a difference by the same mechanism described above for the response to warming. When elevated CO2 is combined with warming, the warming mobilizes soil N, easing N limitation of plant production, and the inhibiting feedback on production associated with higher plant C:N is weakened. This weakening of the C:N feedback allows the direct effects of elevated CO2 to be more strongly manifested, and therefore a stronger response to both elevated CO2 and warming than the sum of the responses to each factor individually. Tissue and Oechell (1987) found the same synergistic effect, resulting in a sustained increase in production with elevated CO2 and warming but only a transient increase with CO2 alone.

Conclusions

Grazing animals can have large effects on ecosystems (Grellman 2002, McLaren and Jefferies 2004, Post and Pedersen 2008, Sjögersten et al. 2008, Rinnan et al. 2009, Cahoon et al. 2011, Kaarlejärvi et al. 2015, Leffler et al. 2019, Min et al. 2021). Our simulations suggest that long‐term exclusion of voles or maintenance of vole populations at high densities can result in large gains or losses of both plants and soil C (Fig. 2). However, the response time of plants and soil to these persistent changes in grazing takes several decades in our model. As a consequence, the full effects of changes in vole densities are never realized when voles cycle on a three‐to‐four‐year time scale. Indeed, cycling at such a high frequency can be incorporated in our model using the long‐term mean density without any substantial change in the predicted long‐term dynamics of the ecosystem. Our simulations excluding and including voles are not purely academic. Recent studies suggest that the Arctic rodent population cycles could dampen in amplitude or be punctuated with periods of non‐cyclic dynamics in response to altered climate conditions, in particular changes in snow conditions (Ims and Henden 2008, Kausrud, 2008, Gilg and Sittler 2009, Brommer et al. 2010, Domine et al. 2018). Our results suggest that the important dynamics for predicting long‐term changes in tundra biogeochemistry in response to climate change are the mean grazer dynamics on decadal scales, not the higher frequency three‐to‐four‐year cycles.

The effects of voles on C and N cycling can have major effects on the biogeochemical responses of tundra to elevated CO2 and warming. These effects need to be explicitly represented in models rather than aggregated with other ecosystem processes. Even if these other processes act in parallel with vole processes, their response to changes in the environment can be very different. Our analysis indicates that failure to explicitly account for voles results in a large underestimation of the responses of tundra to climate warming and to elevated CO2 and warming. Our analysis is analogous to that of Thornton et al. (2007) who found that predicted responses of the terrestrial biosphere to elevated CO2 and climate change was probably overestimated unless N limitation was explicitly represented in models. We recommend that grazing effects be explicitly incorporated when applying models of tundra response to global change.

Our analysis is based on a simple, annual‐time‐step model of ecosystem C and N interactions calibrated to Arctic tundra. The simplicity of the model facilitates causal tracing (Fig. 7) and heuristic analysis of the results, but at the expense of quantitative detail in the dynamics (Rastetter 2017). The results should therefore be confirmed for more complex models with, for example, more detailed representations of vegetation and soil characteristics, finer scale seasonal dynamics, and the effects animals can have on plant‐community composition and soil structure. Although our model was calibrated for Arctic tundra, the qualitative conclusions probably apply more broadly. In our analysis, a key process is vole respiration, which decreases with warming, unlike the increase with warming for plant and microbial respiration. This property is clearly relevant to mammal grazers in cold climates, but not for mammal grazers in warm climates or for insects in any climate. Nevertheless, for these other ecosystems there might be analogous model biases associated with aggregating biogeochemical processes mediated by these grazers with other ecosystem processes. Similarly, the consequences of resource limitation need to be examined for grazers in ecosystems that are bottom‐up regulated. All these possibilities need to be analyzed, first with heuristic models such as the one we use and then incorporated into more detailed biogeochemical models. To perform these analyses, more data such as those in Batzli et al. (1980) and Olofsson et al. (2012) are needed that can be directly applied in these biogeochemical models. Collection of these data will require a biogeochemical, as well as a community, perspective on plant–grazer interactions.

Appendix S1

Click here for additional data file.

Mass‐balance equations
1	dB_Cdt=P_s-L_itC-R_a-G_C	2	dB_Ndt=U_N-L_itN-G_N
3	dD_Cdt=L_itC+L_VC-R_h-Q_CR	4	dD_Ndt=L_itN+U_Nm+L_VN-N_min-Q_NR
		5	dNdt=N_in+N_min+V_UN-U_N-U_Nm-Q_DIN
Allometry and stoichiometry constraints
6	S=B_C∝B_C+1γB_C+1	7	Ψ=B_CB_Nq_B
8	Φ=D_CD_Nq_D
Process equations
Carbon		Nitrogen
9	P_s=g_CSC_aΨk_C+C_aQ10PsT/10	10	U_N=g_NΨSNk_N+NQ10UT/10
11	L_itC=m_CBB_C	12	L_itN=m_NBΨB_N
13	R_a=r_BB_CΨQ10RaT/10	14	U_Nm=g_NmΦD_CNk_Nm+NQ10UmT/10
15	G_C=n_V+g_V-ε_VT-T_0V/10,000	16	G_N=G_C/q_V
17	R_V=r_VG_C-n_VV/10,000	18	V_UN=m_NVV/10,000
19	L_VC=G_C-R_V	20	L_VN=G_N-V_UN
21	R_h=r_DD_CΦQ10RhT/10	22	N_min=m_NmD_NΦQ10mT/10
23	Q_CR=q_DOMQ_NR	24	Q_NR=β_NRD_N
		25	Q_DIN=β_NN