Ecosystems DOI: 10.1007/s10021-016-0034-7  2016 Springer Science+Business Media New York

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations: Implications for a Changing Climate David J. Millar,1,2,3* David J. Cooper,1,2 Kathleen A. Dwire,1,2,4 Robert M. Hubbard,1,2,4 and Joseph von Fischer2,5 1

Department of Forest and Rangeland Stewardship, Colorado State University, Fort Collins, Colorado 80523, USA; 2Graduate Degree Program in Ecology, Colorado State University, Fort Collins, Colorado 80523, USA; 3Present address: Now at: Department of Botany, University of Wyoming, 1000 East University Avenue, Laramie, Wyoming 82071-0333, USA; 4USDA Forest Service Rocky Mountain Research Station, 240 West Prospect Rd, Fort Collins, Colorado 80526, USA; 5Department of Biology, Colorado State University, Fort Collins, Colorado 80523, USA

ABSTRACT Mountain fens found in western North America have sequestered atmospheric carbon dioxide (CO2) for millennia, provide important habitat for wildlife, and serve as refugia for regionally-rare plant species typically found in boreal regions. It is unclear how Rocky Mountain fens are responding to a changing climate. It is possible that fens found at lower elevations may be particularly susceptible to changes because hydrological cycles that control water tables are likely to vary the most. In this study, we fit models of growing season ecosystem-atmosphere CO2 exchange to field-measured data among eight fen plant communities at four mountain fens along a climatic gradient in the Rocky Mountains of Colorado and Wyoming. Differences in growing season net ecosystem production (NEP) among study sites were not well correlated with monsoon precipitation, despite a twofold increase in summer rainfall between two study regions. Our results show

that NEP was higher for fens located at high elevations compared to those found at lower elevations, with growing season estimates ranging from -342 to 256 g CO2-C m-2. This was reflected in the negative correlation of growing season NEP with air temperature, and positive correlation with water table position, as the high elevation sites had the lowest air temperatures and highest water tables due to greater snowpack and later onset of melt. Our results suggest that sustainability of mountain fens occurring at the lower end of the known elevation range may be particularly susceptible to a changing climate, as these peatlands already experience lower snowpack, earlier snow melt, and warmer growing season air temperatures, which are all likely to be exacerbated under a future climate. Key words: peatland; net ecosystem production; mountain; climate change; carbon dioxide; fen.

Received 9 June 2016; accepted 17 July 2016;

INTRODUCTION

Author contributions David Millar: Designed study, performed research, analyzed data, and wrote the paper. David Cooper: Designed study, analyzed data, and wrote the paper. Kathleen Dwire: Designed study, and edited the paper. Robert Hubbard: Designed study, provided data interpretation, and edited the paper. Joseph von Fischer: Designed study *Corresponding author; e-mail: [email protected]

Peatland ecosystems have been sequestering atmospheric carbon (C) for millennia, resulting in an overall global cooling effect due to radiative forcing (Frolking and Roulet 2007). The sequestration of atmospheric C has resulted in the for-

D. J. Millar and others mation of organic soils, the most important defining ecological characteristic of these wetlands, and despite occupying about 3% of the global land surface area, peatland soils now contain approximately half the amount of carbon dioxide (CO2)–C in the atmosphere (Rydin and Jeglum 2006). Carbon dioxide dynamics are the primary control on a peatland’s C balance, and play a key role in their development and persistence. Net ecosystem production (NEP) represents the difference between gross primary production (GPP) and ecosystem respiration (ER). When NEP is positive, it represents a gain of CO2 to an ecosystem, and a loss to the atmosphere when negative (Chapin and others 2006). Peatland CO2 dynamics may be changing, potentially leading to climate feedbacks necessitating a better understanding of the nonlinear dynamics associated with the abiotic and biotic controls on NEP (Belyea 2009). Water table position and air temperature during the growing season act as important controls on CO2 dynamics in peatlands, and are therefore critical driving variables in determining NEP. As water tables drop, more of the soil profile is exposed to aerobic conditions that can increase ER, leading to a reduction in NEP, particularly for groundwaterfed peatlands, known as fens (Chimner and Cooper 2003a; Riutta and others 2007; Wu and others 2013). Recent work has also provided different perspectives; that a decrease in water table depth up to a certain threshold can increase GPP for a Sphagnum-dominated fen (Peichl and others 2014) and a rising water table can reduce NEP by decreasing GPP (Yurova and others 2007). Warmer air temperatures can decrease NEP by increasing ER, particularly when water tables are below the peat surface (Sullivan and others 2007). Further, a rise in air and soil temperature could release ancient, otherwise stable, C from soil organic matter pools through the ‘‘priming’’ of microbes by vegetation in peatlands dominated by grasses, sedges, and shrubs (Walker and others 2016). Fens dominated by wetland plant species, whose main distributions are in boreal regions, exist at their southernmost limit in the Rocky Mountains of western North America (Cooper and Andrus 1994; Chimner and others 2010). These ecosystems provide critical perennially wet habitat for wildlife, support high biodiversity, and sustain rare plant species that are in some cases isolated from their nearest population by over 1000 km (Cooper 1996; Chimner and others 2010). In the Rocky Mountains, temperature decreases, and in most areas, precipitation increases with elevation. Snowmeltderived water provides the hydrologic support for

mountain fens in western North America, thus they rarely occur below the elevation where a large winter snow pack accumulates (Chimner and others 2010). Late-summer precipitation driven by the North American Monsoon is also an important seasonal water source in mountain regions closest to the southwestern US (Costigan and others 2000). Over the past century, no discernible trends in precipitation in the Rocky Mountains have been observed, and models project a range of future trends in annual precipitation (Ray and others 2008). However, the warming climate has led to earlier spring snow melt and runoff, and a reduction in total mountain snowpack in the western US, with more winter precipitation falling as rain instead of snow (Regonda and others 2004; Christensen and Lettenmaier 2007; Ashfaq and others 2013). This shift in the timing of key hydrological processes has reduced late-summer stream flows, leading to deeper water table drawdown, and increasing drought stress in riparian ecosystems (Rood and others 2008). In the southern Rocky Mountains, regionally variable monsoon precipitation further complicates efforts to predict how growing season hydrological cycles are being affected by climate change. It is uncertain whether local and regional hydrological processes can sustain mountain fens, but it is likely that warming air and subsequent changes in watershed hydrological dynamics will affect C cycling in these ecosystems. High elevation mountain ecosystems such as fens in the Rocky Mountains can function as ‘‘sky islands’’, as they are geographically isolated, and surrounded by distinctly different climates at lower elevation that cannot support their existence (Warshall 1994). An upward shift in the elevation range of several species has recently been observed in several mountain regions (Parmesan 2006), particularly at high elevation in western North America (Krajick 2004). A similar shift could also occur for Rocky Mountain fens if warming temperatures coupled with altered hydrological cycles create climatic conditions unsuitable for their sustainability at the low end of their elevation range. Climate change that alters watershed hydrologic regimes and increases air temperature has the potential to reduce annual NEP in mountain fens and affect their long-term sustainability. We hypothesized that such changes are most significant for fens at the low end of their elevation range where they already experience warmer climate conditions and receive less snow. In this study, we developed empirical models of CO2 dynamics driven by meteorological and hydrological variables,

8800 7800 15,300 24,800 10747’45’’W 10758’23’’W 10610’2’’W 10615’15’’W Plant community IDs in parentheses below site names. a Cation and pH estimates are means from two samples, one collected from each plant community, at each site during the growing season of 2013. b Only one measurement was used from a 2006 to 2007 San Juan fen survey (Chimner and others 2010).

2600 3400 2700 3200 5.4b 6.0b 7.0 6.9 SJlow 1 and 2 (Anglica Fen) SJhigh 1 and 2 (Spruce Fen) MBlow 1 and 2 (Sand Lake Fen) MBhigh 1 and 2 (Medicine Bow Peak Fen)

1.7 0.4 20.2 11.6

8.2 2.5 68.8 34.2

2.1 1.6 7.4 0.7

1.9 0.4 3.8 0.3

Elevation (m) Na (ppm) Ca (ppm) Mg (ppm)

Peat chemical propertiesa Site ID

Study Site Locations and Peat Chemical Characteristics Table 1.

We studied eight plant communities at four fens (two communities per fen). Two fens are located in the Medicine Bow Mountains of southern Wyoming, and two in the San Juan Mountains of southwestern Colorado (Table 1). In each study region, we chose one fen located near the lower and one near the upper elevation range of fens (Heidel and Jones 2006; Chimner and others 2010). Mean annual and growing season air temperature among sites decreases with increasing elevation and latitude. Average snowfall increases with elevation and average precipitation during the summer differs between the two study regions. The North American Monsoon has a stronger influence in the San Juan Mountains, where the wettest month of the year is August. Data from two SNOTEL telemetric micrometeorological stations, the Cascade station (ID: 386, elevation: 2710 m) in the San Juan Mountains, and the Brooklyn Lake station (ID: 367, elevation: 3120 m) in the Medicine Bow Mountains were used to provide a long-term (1983–2013) perspective on snow water equivalent (SWE), growing season precipitation, and air temperature in each region. While mosses were present at all sites in varying quantities, the vast majority of fens in the southern Rocky Mountains have Carexderived organic soils (Chimner and others 2002), and our analyses were focused on sedge- and shrub-dominated plant communities. Anglica Fen (SJlow), the low elevation San Juan Mountains study site, occurs in a basin surrounded by granite bedrock outcrops with no surface water inflows or outflows. At this site, we studied one plant community dominated by Carex lasiocarpa (SJlow 1), and another by C. lasiocarpa and Carex utriculata (SJlow 2). Spruce Fen (SJhigh), the high elevation San Juan Mountain site is a sloping fen

K (ppm)

STUDY SITES

pH

Latitude

Longitude

and estimated CO2 fluxes over two growing seasons for a total of eight plant communities at four mountain fens. The fens were located near the low and high ends of their known elevation range in two regions that receive different average amounts of late-summer monsoon precipitation. Our objectives were to: (1) identify trends between growing season NEP in mountain fens and growing season climate and hydrological conditions, (2) use these trends to infer how these ecosystems may be changing in response to climate and (3) compare NEP among mountain fen plant communities occurring at the high and low ends of their elevation range, and in regions with different monsoon influence.

3738’24’’N 3753’53’’N 4120’31’’N 4120’58’’N

Area (m2)

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations

D. J. Millar and others with a surface water inflow and outflow. At this site, we studied one community dominated by Carex saxatilis (SJhigh 1) and another dominated by C. utriculata (SJhigh 2). The low elevation Medicine Bow Mountains study site Sand Lake Fen (MBlow) is gently sloping and surrounded by glacial moraines with no surface water inflow and an ephemeral surface water outlet. We studied a Carex simulata, C. utriculata, and Juncus balticus dominated community (MBlow 1), as well as a Salix planifolia, S. wolfii, and Betula glandulosa community (MBlow 2). The high elevation site was Medicine Bow Peak Fen (MBhigh), a gently sloping fen that receives strong groundwater inflows. At MBhigh, we studied a Carex aquatilis community (MBhigh 1), and a S. planifolia and S. wolfii dominated (MBhigh 2) community.

METHODS Precipitation During the winter months, Moultrie Game Spy cameras were installed at each site and programmed to photograph the majority of the fen, including a centrally located snow depth gage marked at 5-cm increments, to track daily changes in snow depth and determine when each site became snow-free. Peak SWE at each site was determined using snow depth data collected by the onsite cameras. For each site, 2nd order polynomial functions were fit to SWE and snow depth data collected at nearby SNOTEL stations using a Bayesian model and used to model peak SWE as a function of on-site snow depth (Millar 2015). Growing season rain events for each study region were recorded at the low elevation sites using Hobo Onset RG-2 data logging rain gages (Bourne, MA).

CO2 Flux Measurements Fluxes of CO2 were measured in three plots in each plant community. Each plot consisted of an ABS plastic collar, 60 cm 9 60 cm, inserted approximately 5 cm into the soil, and used as the base for a 2.16 9 105 cm3 cubic gas flux chamber. The chamber consisted of an aluminum frame and clear acrylic panels that allowed more than 90% of photosynthetically active radiation (PAR) to pass through. A channel containing closed-cell foam along the bottom of the chamber ensured an airtight seal between the chamber and collars. Because it was not possible to install collars and make chamber-based gas flux measurements in areas dominated by S. planifolia due to the height and

density of the shrubs, for MBlow 2, B. glandulosa was chosen for gas flux plots to represent the shrub stands. Measurements were made in the field using a PP-Systems EGM-4 infrared CO2 gas analyzer (IRGA) (Amesbury, Massachusetts). Concentrations of CO2 were measured at 5-s intervals over a 2-min period and the quadratic slope of the change in concentration over time was used to calculate the flux (Johnson and others 2013). Two battery powered fans within the chamber mixed the air during measurement periods, and on hot midsummer days, ice packs were mounted in the chamber to prevent overheating of chamber air. Measurements of CO2 fluxes were made on an approximately biweekly basis during the growing season of 2011, monthly during 2012, and once in August 2013 at MBlow and SJlow, when the water table was at its lowest position during the study for MBlow. During 2011, mid-day NEP was measured using the clear chamber followed by measurements of ER using a light-proof cover over the chamber. Gross primary production was determined by summing concurrent NEP and ER flux rates (both positive in this case). During the 2012 growing season and in August 2013, fluxes were measured several times during each field day. Measurements using the clear chamber were followed by measures with two shade clothes covering the chamber, reducing PAR by approximately 75 and 50%, followed by an ER measurement using a light-proof cover. This approach allowed for a greater range of PAR used in fitting models to measured fluxes. During each measurement, a sensor attached to the sidewall of the chamber was leveled and used to record PAR and air temperature within the chamber. Water table position was recorded by gaging monitoring wells adjacent to each collar.

Modeling Carbon Dioxide Fluxes We modified equations used by Riutta and others (2007) to model GPP and ER. GPP was modeled for each plant community as a function of PAR and a seasonality term based on a 4-week running average (21 days before, 7 days after) of daily mean air temperature (RAV) [Equation 1]. A rectangular hyperbola function was used to model ecosystem photosynthetic response to incoming PAR, and a Gaussian function was used for the seasonality term, allowing modeled GPP to follow seasonal dynamics associated with plant phenology. h  2 i RAVi RAVoptGPP 0:5 Amax  a  PARi RAVdevGPP GPPi ¼ : ð1Þ e Amax þ a  PARi

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations In Equation 1, Amax (g CO2-C m-2 h-1) represents the asymptotic maximum potential rate of GPP, and a (g CO2-C lmol PAR-1) represents the light use efficiency, or initial slope of the light response function. The parameter RAVoptGPP (C) represents the optimum value of RAV for GPP and RAVdevGPP (C) represents the standard deviation of the Gaussian function, which controls the spread of the distribution. Ecosystem respiration was modeled as a function of air temperature (AT), water table position (WTP), and a seasonality term [Equation 2]. A modified van’t Hoff equation was used to model ER as increasing exponentially with air temperature. The response of ER to water table position was modeled as a negative exponential equation, and a Gaussian function similar to that of the GPP model was used to account for seasonal variation in ER. h  2 i   ATi 10 10

ERi ¼ R10  Q10

 ebWTPi  e

0:5

RAVi RAVoptER RAVdevER

: ð2Þ

In Equation 2, R10 (g CO2-C m-2 h-1) represents ER at 10C when other model factors are not limiting, Q10 represents the rate of increase in ER per 10C increase in air temperature, b (g m-2 cm-1) represents the initial slope of the rate of increase in ER per decrease in water table position below the peat surface. RAVoptER (C) and RAVdevER (C) represent the optimum RAV value for ER and the standard deviation of the Gaussian function controlling seasonality in ER, respectively.

Model Development and Evaluation We fit models to the measured data using Bayesian methods in R statistical software (R Development Core Team, Austria). Model parameters were estimated using Markov chain Monte Carlo (MCMC) analysis in the rjags package for R (Plummer 2011). A total of 100,000 iterations were used with 4 MCMC chains, after a burn-in of 30,000 iterations. Vague priors were used for all model parameters. Uniform distributions with limits ranging between 0 and 30 were used for RAVoptGPP, RAVdevGPP, RAVoptER, and RAVdevER priors, uniform distributions with limits ranging between 0 and 0.1 were used for a priors, gamma distributions with shape and rate parameters equal to 0.001 were used for Amax, rprocGPP, R10, Q10, and rprocER priors, and a beta distribution with both shape parameters equal to 1 was used for the b priors. We used meanweighted variance in modeling GPP and ER to account for variance that increased with both flux rates (Equations 3 and 4).

 NormalðlGPPi ; rprocGPP  lGPPi Þ GPPobs i

ð3Þ

ERobs  NormalðlERi ; rprocER  lERi Þ: i

ð4Þ

Equations 3 and 4 represent the likelihood functions for the observed GPP (GPPobs i ) and ER obs (ERi )measured in each plant community. In the equations, lGPPi and lERi represent the predicted values of GPP and ER, respectively, and rprocGPP and rprocER represent the process variance associated with those predictions.

Modeling Growing Season CO2 Fluxes Once model parameters were estimated, GPP and ER models were run for all sites, for the period May 28 through September 19 in 2012 and 2013. Models were run using hourly PAR, air temperature, and water table position measured throughout both growing seasons. Continuous PAR measurements were recorded using a Campbell Scientific CR10X data logger (Logan, UT) equipped with an Apogee Instruments SQ-110 quantum sensor (Logan, UT). PAR measurements used to drive growing season GPP models for the San Juan sites and Medicine Bow sites were made at the respective low elevation sites in each of these regions. At each site, an In-Situ Barotroll logger (Fort Collins, CO) was used to record hourly air temperature. In-Situ Rugged Troll pressure transducers were used to measure hourly water table position in monitoring wells installed in each plant community. Because of the patchy distribution of shrub and graminoid communities at MBlow, and no obvious hummock-hollow microtopography, one pressure transducer was used for both communities. During the study, the water table fell below the elevation of the monitoring well pressure transducer at SJlow 2 and SJhigh 2. This resulted in several periods without data, and models were run with the water table at -54 cm from 7/20/12 to 7/ 31/12, 8/20/12 to 9/18/12, 6/28/13 to 6/29/13, and 7/12/13 to 9/18/13 for SJlow 2, and -16 cm from 6/ 24/13 to 7/3/13, 8/11/13 to 8/22/13, 9/4/13 to 9/ 11/13, and 9/15/13 to 9/18/13 for SJhigh 2. Hourly and cumulative growing season GPP and ER were estimated using MCMC, with a total of 6000 iterations after a burn-in of 3000 iterations. Hourly and cumulative growing season NEP estimates were determined similarly, as the difference between GPP and ER. Differences in simulated cumulative growing season NEP were compared between all sites using robust Bayesian Estimation,

D. J. Millar and others using 1000 randomly selected MCMC samples from each posterior distribution being compared (Kruschke 2013). Once growing season CO2 flux estimates were determined, we developed a Bayesian Multiple Linear Regression (BMLR) model using mean growing season air temperature and water table position as predictor variables and mean growing season NEP as the response variable (Equation 5). PAR is also an important driver of NEP, and increased cloudiness has been shown to significantly decrease ecosystem C uptake in a subarctic fen by reducing incoming PAR by approximately 50% (Nijp and others 2015). However, despite having less precipitation than the San Juan Mountains, ecosystems in the Medicine Bow Mountains typically experience similar cloudiness accompanied by a reduction in incoming PAR late in the afternoon during the growing season (Sanchez and others 2015). We compared mean daytime PAR (08:00–18:00) between the two regions using robust Bayesian estimation (Kruschke 2013), and found that it varied by only 3% (difference in mean growing season PAR = 34.7 lmol m-2 s-1, 95% CI = 66.9, 0.33). Therefore, we did not include mean growing season PAR in the BLMR model. NEPj ¼ b0 þ b1 ATj þ b2 WTPj :

ð5Þ

NEP is the cumulative growing season NEP, determined by summing hourly flux estimates, AT is mean growing season air temperature (C), and WTP is mean growing season water table position (cm). b0, b1, and b2 represent regression coefficients, all of which had normal distributions for vague priors with shape parameters of 0 and 1 9 10-6 for mean and variance, respectively. Equation 6 represents the likelihood function for growing season NEP (NEPmodeled ) for each plant j community, derived using the previously described deterministic models. NEPmodeled  NormalðlNEPj ; rprocNEP Þ: j

ð6Þ

In Equation 6, lNEPj represents the predicted values of growing season NEP using the BMLR model, and rprocNEP represents the process variance associated with those predictions. Parameters for the BMLR models were also estimated using MCMC, with a total of 25,000 iterations with 4 chains, after a burn-in of 10,000 iterations. We used the Gelman-Rubin diagnostic to test model convergence for all parameters (Gelman and Rubin 1992).

RESULTS Growing Season Climate and Hydrological Conditions Mean growing season air temperature decreased with increasing elevation across sites during both study years (Table 2). The warmest mean growing season air temperature, 15.4C, occurred at SJlow in 2012, while the coldest occurred at the highest site, SJhigh, 10.5C in 2012 and 2013. With the exception of MBhigh in 2012, the high elevation sites had greater peak SWE than either of the low elevation sites. The low peak SWE estimated for MBhigh in 2012 was likely a function of wind redistribution of snow, as this site was much more exposed to high winds than the other three, and peak SWE at the nearby SNOTEL site (North French Creek, #668) at an elevation of 3100 m, was 48.5 cm. During both study years, the San Juan Mountain sites had higher peak SWE than their high and low elevation counterparts in the Medicine Bow Mountains. Within each region, the low elevation sites were snow-free from several days to approximately a month earlier than the high elevation sites (Table 2). In 2012, both Medicine Bow sites were snow-free by April 2, more than a month earlier than in 2013, and earlier than both San Juan sites in 2012. Peak SWE in both 2012 and 2013 was below the 20-year average in both regions. In particular, 2012 had the 3rd lowest peak SWE since 1983 in the Medicine Bow Mountains. June through August precipitation was below average in both regions during the study years; however deviation from the mean was most significant in the Medicine Bow Mountains where rainfall was approximately half of the 20-year average. Both SNOTEL sites showed a trend in increasing mean annual air temperature of approximately 0.2C y-1, from 1990 to 2013 (0.62 £ R2 £ 0.67). Thus, mean growing season air temperatures among all sites were higher than average for the prior two decades. The San Juan Mountains received 77–94% more summer (June–August) precipitation than the Medicine Bow Mountains. June precipitation ranged from 5 to 16 mm for both regions. However, late-summer precipitation (July–August) in the San Juan Mountains was 40–71 mm in 2012 and 2013, compared with 17–40 mm for the Medicine Bow Mountains. The 2013 growing season was wetter than 2012 in both regions, in part due to unusually high precipitation during September (Table 2; Figure 1).

-25.4 (47.0), 17.7 (82.5) -43.0 (44.1), 31.6 (78.7) 188.8 (74.5), 234.4 (58.3) 188.2 (67.6), 256.2 (51.9) 210.0 (66.2), 175.6 (54.5) -341.8 (167.2), -53.4 (127.8) 70.7 (119.9), 237.8 (55.6) 150.5 (977.4), 234.6 (53.6) -13.9 (8.8), -40.0 (16.5) -15.3 (6.9), -31.4 (16.6) -4.2 (12.7), -12.6 (3.4) -3.8 (2.2), -6.2 (0.3) -17.6 (9.7) -52.0 (30.4) -24.1 (16.5), -12.4 (4.8) -5.9 (3.7), -4.8 (1.3)

Growing season water tables at the high elevation plant communities averaged 19 cm closer to the soil surface (mean depth range: 3.8–24.1 cm) than the low elevation communities (mean depth range: 13.9–52.0 cm) (Table 2). Within both regions, low elevation sites had the earliest and most dramatic water table decline (Figure 1). Rain events, which began at the end of June, appeared to play a more significant role in maintaining near surface water tables in the low elevation plant communities in both regions, by periodically reversing steep declines. Both high elevation communities produced less ‘‘flashy’’ hydrographs, and were in general more temporally stable than those at the low elevations.

Gross Primary Production and Ecosystem Respiration

Values are means with standard deviation of the posterior estimates in parentheses. Values in parentheses represent standard deviation. c Values in parentheses represent standard deviation of the posterior estimate. b

a

MBhigh

MBlow

SJhigh

Mar-4 Mar-10 Mar-21 Apr-17 Feb-23 Apr-17 Feb-14 Apr-17 2012 2013 2012 2013 2012 2013 2012 2013 SJlow

32 30 48 46 34 34 24 48

(1.5) (1.0) (0.74) (0.84) (0.37) (0.37) (0.38) (0.30)

Apr-5 Apr-21 May-15 May-22 Mar-31 May-15 Apr-2 Jun-1

124 224 124 224 63 199 63 199

15.4 14.9 10.5 10.5 13.0 12.5 11.2 10.7

(6.7) (6.3) (5.3) (5.4) (7.1) (6.9) (5.8) (5.5)

Mean growing season NEP (g CO2-C m-2 gs-1)c 1, 2 Mean growing season water table (cm)b 1, 2 Mean growing season air temperature (C)b Growing season precipitation (mm) Snow-free date Peak SWEa (cm) Date of peak SWE Year Site

Table 2. Climate and Hydrological Characteristics and Growing Season Net Ecosystem Production at Each Site

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations

Modeled GPP and ER generally matched measured fluxes well, however three of the four ER models for the high elevation plant communities had low R2 values, around 0.30 (Table 3; Figure 2). The range of air temperatures recorded during CO2 flux measurements at the high elevation sites were, on average, 11C lower than the range at the low elevation sites, and the range in recorded water table positions were, on average, also lower (11 cm), which may explain the low R2 values for the ER models at these sites. Nevertheless, all predicted CO2 fluxes showed homoscedastic variance across plant communities (Figure 2). Further, measured GPP and ER fluxes showed seasonal trends consistent with the seasonality terms in each model (Figure 3). Field-measured GPP rates ranged from 0.003 to 1.26 g CO2-C m-2 h-1 across measurement years, and increased with elevation within each region and with latitude between the high and low elevation communities in each region (Figure 2). Similarly, cumulative growing season GPP generally increased with elevation and latitude (Figure 4). During both study years, cumulative growing season GPP was lowest for SJlow 1, 373 g CO2-C m-2 in 2012 and 350 g CO2-C m-2 in 2013. Highest cumulative growing season GPP estimates during both study years were observed in MBhigh 2, ranging from 780 g CO2-C m-2 -in 2012 and 741 g CO2-C m-2 in 2013. Parameter estimates from the GPP models show similar light use efficiency and seasonality among plant communities, with the exception of SJlow 1 and SJlow 2, which had much lower estimates of Amax, both approximately half that of the other plant communities (Table 3).

D. J. Millar and others Fig. 1. Mean daily growing season water table dynamics from centrally located monitoring wells for high and low elevation sites in the San Juan and Medicine Bow Mountains. Daily precipitation events for the San Juan (black) and Medicine Bow (gray) Mountains are represented as vertical bars.

Field-measured ER rates ranged from 0 to 0.53 g CO2-C m-2 h-1 across measurement years, and followed a similar trend as measured GPP, increasing with elevation within regions and with latitude between high and low elevation plant communities between regions (Figure 2). However, unlike cumulative growing season estimates of GPP, cumulative growing season ER was more variable among communities, both within and among sites (Figure 3). The lowest cumulative growing season ER occurred in the two SJhigh communities during both years, ranging from 304 to 348 g CO2-C m-2 for both plant communities. Highest cumulative growing season ER occurred in both years at the Medicine Bow sites. In 2012, the highest ER occurred at MBhigh, 596 g CO2-C m-2 for MBhigh 1 and 542 g CO2-C m-2 for MBhigh 2. In 2013, highest cumulative growing season ER occurred at MBlow, 960 g CO2-C m-2 for MBlow 1 and 629 g CO2-C m-2 for MBlow 2. Parameter estimates for the ER models revealed a similar response in ER to air temperature across plant communities, but a varied response in ER to changes in water table position (Table 3).

Net Ecosystem Production With the exception of the plant communities at MBlow in 2012, the probability of plant communities acting as C sinks during the growing season was higher for the high elevation plant communities than that of the low elevation ones (Figure 4). Cumulative growing season NEP was similar be-

tween years for the San Juan plant communities at both elevations. In the Medicine Bow Mountains, MBlow communities had much lower NEP in 2013 than 2012, driven by a lower water table in 2013. Conversely, cumulative growing season NEP at MBhigh 2 was similar between years, whereas MBhigh 1 increased approximately twofold between 2012 and 2013, which corresponded with an increase in mean growing season water table position between years (Table 2; Figure 3). The lowest cumulative growing season NEP in 2012 occurred at SJlow, with SJlow 1 representing the only negative NEP estimate, -25 g CO2-C m-2, for this during this study year. SJlow 2 represented the lowest positive cumulative growing season NEP at 18 g CO2-C m-2. In 2013, the lowest cumulative growing season NEP occurred at the two MBlow plant communities, -342 g CO2-C m-2 for MBlow 1 and -53 g CO2-C m-2 for MBlow 2. In addition, NEP was lower for both plant communities at SJlow in 2013, with rates of -43 g CO2-C m-2 for SJlow 1 and 32 g CO2-C m-2 for SJlow 2. Highest growing season NEP occurred at the high elevation sites in both years of the study, 237 g CO2-C m-2 for MBhigh 2 in 2012, and 256 g CO2-C m-2 for SJhigh 2 in 2013 (Table 2; Figure 3). Results from the Bayesian estimation analysis revealed that differences between mean cumulative growing season NEP were credible for the vast majority of community-to-community comparisons, with the exception of MBhigh 2 and SJhigh 2 in 2012, and MBlow 2 and SJlow 1 in 2013 (Table 4). Cumulative growing season NEP for the high ele-

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations Table 3.

Gross Primary Production and Ecosystem Respiration Model Parameter Estimatesa SJlow 1

SJlow 2

SJhigh 1

SJhigh 2

Amax a RAVoptGPP RAVdevGPP rGPP

0.57 (0.13) 0.0018 (0.00097) 23.8 (3.9) 10.5 (2.4) 0.40 (0.038)

0.57 (0.13) 0.011 (0.0062) 24.4 (3.8) 12.0 (3.3) 0.67 (0.066)

1.40 (0.46) 0.0030 (0.0012) 23.2 (4.5) 10.6 (2.2) 0.52 (0.047)

1.24 (0.59) 0.0070 (0.0054) 22.8 (4.9) 10.3 (3.0) 0.85 (0.014)

N R2

78 0.65

69 0.49

73 0.72

69 0.61

Q10 R10 b RAVoptER RAVdevER rER

1.2 (0.088) 0.14 (0.032) 0.0069 (0.0037) 21.4 (4.4) 10.3 (3.6) 0.12 (0.13)

1.3 (0.099) 0.15 (0.035) 0.003 (0.0023) 19.7 (4.7) 10.8 (5.6) 0.22 (0.025)

1.5 (0.30) 0.14 (0.051) 0.016 (0.01) 19.8 (6.5) 17.7 (6.2) 0.25 (0.027)

1.2 (0.14) 0.15 (0.051) 0.012 (0.0062) 21.0 (5.9) 15.4 (5.8) 0.13 (0.014)

N R2

58 0.74

53 0.75

49 0.32

52 0.30

MBlow 1

MBlow 2

MBhigh 1

MBhigh 2

Amax a RAVoptGPP RAVdevGPP rGPP

1.32 (0.29) 0.0021 (0.00057) 22.9 (4.1) 10.8 (2.0) 0.68 (0.056)

1.24 (0.25) 0.0022 (0.00062) 24.0 (3.7) 11.3 (1.9) 0.66 (0.053)

1.35 (0.38) 0.0027 (0.00084) 23.2 (4.5) 12.3 (2.6) 0.86 (0.077)

1.59 (0.49) 0.0035 (0.0012) 21.9 (4.9) 10.3 (2.4) 1.1 (0.09)

N R2

97 0.80

101 0.80

85 0.72

88 0.75

Q10 R10 b RAVoptER RAVdevER rER

1.5 (0.2) 0.12 (0.032) 0.019 (0.0024) 24.4 (3.9) 18.0 (3.9) 0.15 (0.014)

1.3 (0.096) 0.15 (0.030) 0.0095 (0.0028) 24.2 (4.0) 18.8 (4.1) 0.16 (0.014)

1.2 (0.15) 0.17 (0.038) 0.012 (0.0058) 18.7 (6.7) 19.7 (5.9) 0.25 (0.027)

1.2 (0.075) 0.19 (0.032) 0.04 (0.0072) 21.6 (5.2) 17.5 (5.1) 0.15 (0.016)

N R2

71 0.67

75 0.58

52 0.30

51 0.75

a

Standard deviation of the posterior parameter estimates are shown in parentheses.

Fig. 2. Observed versus posterior mean predicted ecosystem respiration and gross primary production for all study sites.

D. J. Millar and others

Fig. 3. Seasonal trends of measured ecosystem respiration and gross primary production fluxes for all years of the study.

vation plant communities was substantially higher than that of the low elevation communities (>150 g CO2-C m-2 gs-1) for both study years in the San Juan Mountains, and in 2013 for the Medicine Bow Mountains communities. This trend was not the case for the Medicine Bow Mountain communities in 2012, due in part to the particularly low cumulative growing season NEP at MBhigh 1, and high cumulative growing season NEP for both communities at MBlow. In general, comparisons between communities from SJhigh and MBhigh revealed the lowest differences in cumulative growing season NEP, with the exception of MBhigh 1 versus SJhigh 1 and 2 in 2012. Alternatively, MBlow plant communities had substantially higher cumulative growing season NEP than the communities of SJlow in 2012, while this trend was somewhat reversed in 2013.

Fig. 4. Mean cumulative growing season estimates of ecosystem respiration (gray bars), gross primary production (white bars), and net ecosystem production (black bars) for eight plant communities during each study year. Error bars represent standard deviation of the posterior estimate. Dashed lines divide the San Juan sites from the Medicine Bow sites, and the dotted lines separate high and low elevation sites within each region.

Fig. 5. Posterior distributions for cumulative growing season net ecosystem production (NEP) for each plant community (two per fen). Negative NEP represents a loss of carbon and positive NEP represents a gain of carbon.

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations Table 4. Differences Between Mean Cumulative Growing Season Net Ecosystem Production Among Plant Communities During Both Study Years 2012

2013

Plant communities being compared

Difference in means (95% HDI)a (g CO2-C m-2 gs-1)

SJlow 1 SJlow 1 SJlow 1 SJlow 2 SJlow 2 SJlow 1 SJlow 2 SJlow 1 SJlow 2 MBhigh 1 SJlow 2 MBhigh 1 MBhigh 1 MBhigh 1 SJlow 1 MBhigh 1 SJlow 2 MBlow 2 MBlow 2 MBlow 2 SJhigh 1 SJhigh 1 SJhigh 1 SJlow 1 MBlow 1 MBlow 2 MBlow 1 SJhigh 2

-270.0 -265.1 -253.6 -231.2 -226.9 -218.7 -217.6 -205.8 -176.4 -164.3 -163.0 -160.0 -151.8 -126.4 -100.2 -92.5 -68.8 -63.1 -59.0 -52.7 -45.9 -39.9 -36.0 -35.8 -15.0 -12.2 -10.9 -2.8

MBhigh 2 SJhigh 2 MBlow 1 SJhigh 2 MBhigh 2 SJhigh 1 MBlow 1 MBlow 2 SJhigh 1 SJhigh 2 MBlow 2 MBhigh 2 MBlow 1 SJhigh 1 MBhigh 1 MBlow 2 MBhigh 1 MBhigh 2 SJhigh 2 MBlow 1 MBhigh 2 SJhigh 2 MBlow 1 SJlow 2 MBhigh 2 SJhigh 1 SJhigh 2 MBhigh 2

(-273.0, -265.4) (-269.7, -260.7) (-258.0, -248.2) (-237.2, -224.7) (-234.0, -220.9) (-224.1, -213.5) (-223.9, -210.5) (-210.2, -201.2) (-183.5, -169.3) (-172.3, -156.2) (-169.3, -157.0) (-167.6, -151.0) (-160.5, -144.1) (-135.1, -117.9) (-108.4, -92.8) (-99.4, -83.8) (-77.5, -59.8) (-68.0, -58.4) (-63.9, -54.1) (-57.9, -47.2) (-51.6, -40.1) (-46.1, -34.4) (-41.5, -28.7) (-42.3, -30.0) (-20.6, -9.4) (-17.6, -6.7) (-16.0, -5.0) (-7.8, 2.0)b

Plant communities being compared

Difference in means (95% HDI) (g CO2-C m-2 gs-1)

MBlow 1 MBlow 1 MBlow 1 MBlow 1 MBlow 1 MBlow 2 SJlow 1 MBlow 1 MBlow 1 MBlow 2 SJlow 1 SJlow 2 MBlow 2 SJlow 1 SJlow 2 MBlow 2 SJlow 1 SJlow 2 SJlow 2 MBhigh 1 MBhigh 1 MBlow 2 SJhigh 1 SJlow 1 SJhigh 1 MBhigh 1 MBhigh 2 MBlow 2

-600.6 -570.6 -543.3 -494.1 -354.9 -305.8 -300.3 -300.3 -297.6 -281.2 -278.5 -237.5 -236.4 -234.0 -217.9 -199.3 -198.5 -169.8 -139.8 -102.4 -80.5 -70.7 -68.0 -62.5 -48.9 -33.0 -23.0 -3.6

SJhigh 2 MBhigh 2 SJhigh 1 MBhigh 1 SJlow 2 SJhigh 2 SJhigh 2 SJlow 1 MBlow 2 MBhigh 2 MBhigh 2 SJhigh 2 SJhigh 1 SJhigh 1 MBhigh 2 MBhigh 1 MBhigh 1 SJhigh 1 MBhigh 1 SJhigh 2 MBhigh 2 SJlow 2 SJhigh 2 SJlow 2 MBhigh 2 SJhigh 1 SJhigh 2 SJlow 1

(-611.5, -589.7) (-581.5, -559.7) (-554.2, -529.6) (-505.1, -483.2) (-365.8, -344.0) (-314.0, -297.6) (-303.0, -294.8) (-311.2, -289.4) (-308.5, -283.9) (-289.4, -272.7) (-281.2, -273.0) (-243.5, -231.8) (-245.2, -226.6) (-239.4, -229.3) (-223.9, -211.8) (-208.8, -190.3) (-204.2, -193.0) (-175.8, -163.0) (-146.1, -131.6) (-108.4, -96.4) (-86.3, -74.3) (-80.8, -60.9) (-74.0, -63.1) (-68.5, -57.3) (-54.3, -43.4) (-39.9, -26.8) (-28.4, -18.9) (-11.2, 5.6)b

a

Lower and upper limits of the Highest Density Interval (HDI) are shown in parentheses Denotes a difference in mean cumulative growing season NEP that is not credible

b

The BMLR model revealed that cumulative growing season NEP decreased with increasing mean growing season air temperature, and with decreasing mean growing season water table position (R2 = 0.67) (Figure 6, Table 5). Cumulative growing season NEP among plant communities generally increased with elevation (Figures 3, 5). This was particularly true for the San Juan fens where, on average, NEP was an order of magnitude higher at SJhigh than SJlow in both years. NEP was an order of magnitude higher at MBhigh than MBlow in the Medicine Bow in 2013, while in 2012, NEP was slightly higher at MBlow than MBhigh.

DISCUSSION Net Ecosystem Production and CO2 Dynamics Among Study Fens Our results show a trend of growing season NEP increasing with elevation among the eight plant communities studied. Our BLMR results show this was driven by a strong positive correlation with water table elevation and negative correlation with air temperature, both key abiotic factors in controlling peatland CO2 dynamics (Sullivan and others 2007). Growing season NEP ranged from negative to positive among plant communities,

D. J. Millar and others

Fig. 6. Mean growing season net ecosystem production (NEP) versus mean growing season air temperature (A) and mean growing season water table position (below peat surface) (B). Error bars represent standard deviation of the posterior estimates of NEP.

Table 5. Mean Parameter Estimates for the Bayesian Multiple Linear Regression Model Parameter estimatea b0 b1 b2 rprocNEP

438.7 -16.7 7.0 103.9

(175.8) (15.2) (1.9) (22.2)

a Values in parentheses represent standard deviation of the posterior parameter estimate.

across a temperature range of a few degrees C, which is consistent with recent findings from a subarctic fen in northern Sweden (Wu and others 2013). Likewise, the decrease in NEP with declining water table position observed in this study is similar to that in boreal peatlands of Europe and North America (Bubier and others 2003; Riutta and others 2007), and other fens in the Rocky Mountains (Chimner and Cooper 2003b; Schimelpfenig and others 2014). Further, posterior estimates of cumulative growing season NEP with standard deviations overlapping zero, and negative means were found at both of the low elevation sites during the study period. This supports our hypothesis that

continued increases in air temperature and the subsequent impacts on water table dynamics have the potential to threaten the long-term sustainability of mountain fens found at low elevations. Our analyses indicate that the low elevation fen plant communities can function as CO2 sources, or only marginal CO2 sinks. The high elevation fens functioned as CO2 sinks, with a mean cumulative growing season NEP of 195 g CO2-C m-2 for the four plant communities analyzed. Cumulative growing season NEP for the SJhigh and MBhigh plant communities was similar to those reported for higher latitude peatlands with similar vegetation (Riutta and others 2007; Adkinson and others 2011; Maanavilja and others 2011) and in some cases higher (Nilsson and others 2008). The low elevation fens generally had much lower growing season NEP, with net losses of CO2 for one or both of the study years, and an overall mean NEP rate of about 1 g CO2-C m-2, similar to CO2 fluxes reported for hydrologically modified fens in the Rocky Mountains (Chimner and Cooper 2003b). Although many of the GPP model parameters were similar across sites, the substantially lower estimates of Amax- at SJlow may be related to sparse cover of aboveground vegetation, with both communities being dominated by C. lasiocarpa, compared to higher vegetative cover in the other plant communities. These differences limited growing season GPP at this site, which contributed to low overall growing season NEP estimates that had standard deviations overlapping zero for both communities during both growing seasons. Such limitations associated with photosynthetic uptake, in addition to the different sensitivities in fen ER to water table decline, further complicated growing season NEP comparisons related to monsoon precipitation among mountain fen plant communities. The estimates for Q10 among sites were similar; however these estimates were lower than those typically cited in the literature for ER models (Davidson and others 2006). Unlike many other temperature-dependent models of ecosystem and soil respiration (Davidson and others 1998; Melillo and others 2011; Cable and others 2013; Ricker and others 2014), our model used air temperature rather than soil temperature as a driving variable. Because air temperature varies more on a diurnal basis than soil temperature, our ER models were less sensitive to temperature than others have reported. Low R2 values for observed ER fluxes versus posterior means of the predicted fluxes at three of the four high elevation plant communities were lower than those of the other plant communities in

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations this study. These results are similar to those found by Schneider and others 2012 for a minerogenic (groundwater fed) hollow in a northeastern Europe peatland complex, where a water table close to the peat surface limited ER fluxes, and led to low R2 values for an ER model driven by air temperature and water table position. Further, the uncertainty associated with posterior distributions of simulated ER fluxes was similar across sites. Our study included eight plant communities distributed among four mountain fens, at two elevations within the San Juan Mountains of southwest Colorado and the Medicine Bow Mountains of southern Wyoming. However, the dominant plant species among these communities are common in Rocky Mountain fens, and representative of fen ecosystems across the region (Heidel and Jones 2006; Chimner and others 2010). Further, in most cases, GPP and ER model parameters were similar among sites, suggesting that differences in species composition did not necessarily translate to differing responses in ecosystem-scale CO2 fluxes. While our analyses spanned a range in elevation of 800 m (from 2600 to 3400 m) and a range in latitude of approximately 400 km, our results remain limited to the climatic and hydrogeological conditions present at the study sites.

Annual CO2 Fluxes The growing season represents a relatively small part of the year in montane and subalpine environments. Wintertime fluxes of CO2 have been identified as an important contribution to the overall annual CO2 balance of boreal peatlands (Aurela 2002), and nearly half of the C gained through NEP can be lost through ER outside of the growing season Sagerfors and others (2008). Mast and others (1998) reported an average wintertime NEP rate of -7.2 9 10-3 g CO2-C m-2 h-1 for a subalpine fen dominated by C. aquatilis and C. utriculata, in the Colorado Rocky Mountains. Applying this average flux rate to the remaining 250 days of the year outside the growing season study period results in a total NEP flux of 49 g CO-2-C m-2 winter-1. NEP on an annual basis is lower than growing season NEP, and including winter fluxes has critical implications for low elevation fens that have negative growing season NEP, as we found in this study. In particular, annual NEP for the two SJlow plant communities would be negative, acting as a CO2 source during both years of this study. Fens in the southern Rocky Mountains under natural conditions have been reported to act as

both sinks and sources of CO2 on an annual basis (Wickland and others 2001; Chimner and Cooper 2003b), and variability in temperature and precipitation can directly and indirectly affect their source/sink role. Although the results of this study represent only two growing seasons, in any year certain patterns will persist: in general, air temperature decreases with elevation and snowpack increases with elevation, as is the case in our study areas. Thus, the elevation-related trends in NEP for the mountain fen plant communities will likely be true for future years, regardless of year-to-year variability in temperature and precipitation.

Elevation and Monsoon Effects on Fen Net Ecosystem Production We were able to capture trends in peak SWE and the timing of snow melt associated with the elevations that support fens in the western US (Moore and others 2014). The low elevation fens had lower peak SWE and earlier snowmelt due to warmer early summer air temperatures, resulting in a longer snow-free season than the higher elevation sites. Earlier snow melt led to an earlier growing season water table decline at both low elevation fens, even though they exist in different hydrogeological settings. This contributed to higher ER relative to GPP, and lower NEP compared to the high elevation sites that maintained shallow water tables throughout most of the growing season. We estimated growing season CO2 fluxes from the end of May through mid-September during both study years at all sites, regardless of when they became snow-free. This allowed the comparison of CO2 fluxes among sites during the same time period in both years, but did not account for the period between melt out and the start of our model simulations, which differed by over a month in some cases between high and low elevation sites. However, NEP in the period between snow melt and the start of the model simulations was likely low or possibly negative due to low plant biomass in the spring (Blanken 2014), as well as lower daily PAR and higher rates of nighttime ER associated with shorter day length earlier in the year (Wohlfahrt and others 2013). Although the San Juan sites received almost double the growing season rainfall as the Medicine Bow sites, there were no discernible differences in fen NEP between these two regions that were directly associated with summer precipitation. While mean growing season water table position was an important predictor of cumulative growing season NEP among the plant communities in this study,

D. J. Millar and others the response of water tables to precipitation events can be spatially and temporally variable in wetlands due to complex factors such as local topography and geomorphology (Tufford 2011; Vidon 2012). Cumulative growing season NEP estimates at the SJhigh plant communities were strongly positive due in part to relatively high, stable water tables. The water table in these communities varied little during the study years and did not respond strongly to precipitation events. At SJlow, which received the same amount of growing season precipitation as SJhigh, cumulative growing season NEP estimates were much lower, due in part to higher ER resulting from deeper water tables. Unlike SJhigh, the water table rose rapidly followed precipitation events at SJlow, where rain followed dry periods when the water table steadily declined. Despite receiving considerably less rain from July through August, the water table responded similarly in Medicine Bow plant communities, with a rapid rise following rain events at MBlow, and subtle water table rises at MBhigh. The disparity in the responses of water tables to precipitation between the high and low elevation sites is likely due to the position of the water table relative to the peat surface, prior to rain events. Saturated hydraulic conductivity (Ks) decreases with depth in peat soils, with surface peat layers being highly conductive, relative to the more decomposed older peat layers lower in the soil profile (Letts and others 2000; Schimelpfenig and others 2014). For plant communities in sloping fens with shallow water tables, like SJhigh and MBhigh, it is possible for infiltrated precipitation to move laterally through surface peat layers relatively quickly, dampening the magnitude of water table rise. Further, soil porosity decreases with depth in peat (Letts and others 2000; Schimelpfenig and others 2014), increasing the magnitude of water table rise to infiltrated precipitation at lower depths. Although this was the case, since water table decline occurred earlier in the low elevation sites and the water table was considerably lower than that of the high elevation sites during the monsoon season, rises in water table were short-lived, and insufficient to maintain saturated soil conditions near the peat surface, comparable to the high elevation sites. Despite below average summer precipitation during both study years at MBlow the mean growing season water table was 18 cm below the soil surface in 2012, but substantially lower in 2013, at -52 cm. The relatively high mean growing season water table position in 2012 despite low summer precipitation rates may have been due to higher than average groundwater flow into the fen. In the

winter of 2010/2011, above average snowfall in the Medicine Bow Mountains contributed 1.3 m of precipitation during the 2011 water year, which is approximately 44% above the long-term average annual precipitation. Because MBlow is fed primarily by groundwater, with no surface water inflows, it is possible that the effects of above average groundwater recharge during the 2011 water year persisted through the 2012 growing season, causing the water table to remain closer to the surface despite sparse rainfall. By the 2013 growing season, the water table was dramatically lower, and in a similar manner as 2012, may result from a legacy effect of groundwater flow from 2 years prior. The importance of time lags in understanding downgradient responses to upgradient groundwater recharge has been recently noted (Tesoriero and others 2015; Van Meter and Basu 2015). Our results suggest that time lags should be considered when determining the influence of snowmelt on mountain fen hydrological dynamics, particularly for fens lacking surface water inflow.

Rocky Mountain Fens in a Future Climate The strong climate gradient between low and high elevation fen plant communities in the southern Rocky Mountains controlled differences in cumulative growing season NEP, as shown in the Bayesian estimation and BMLR model results, regardless of plant community composition. In most cases, the largest difference in cumulative growing season NEP between plant communities was between high and low elevation sites. Over the coming decades, average annual temperatures are expected to increase by 1.0–2.4C in the southern Rocky Mountains (Christensen and others 2004) and reductions in snowpack of 10–40% have been predicted for elevations between 2500 and 3000 m in Colorado (Christensen and Lettenmaier 2007) an elevation range that includes the lower elevation limit for fens in the region (Chimner and others 2010). The changing climate is therefore likely to alter the hydrological regime of mountain fens during the growing season, lowering water tables (Rood and others 2008) and consequently decreasing NEP, particularly for low elevation fens. We observed the lowest cumulative growing season NEP in fen communities at the low end of their known elevation range in both mountain regions. Reductions in NEP associated with a warming climate may convert low elevation fens, with already low NEP, from net sinks to net sources of CO2. Therefore, the long-term sustainability of fens at low to middle elevations in the Rocky

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations Mountains, and other mountain ranges in the western U.S. may be in jeopardy, where their organic soil could be lost through decreased NEP. Furthermore, mountain fens containing sedge-derived peat, like the fens in this study, may be particularly susceptible to increased decomposition resulting from lowered water tables, since most of the CO2 fixed through GPP is allocated belowground (Chimner and others 2002). An important consideration in understanding the long-term sustainability of mountain fens to climate change is their hydrogeological and topographic setting. Such factors may have contributed to the differences in cumulative growing season NEP between the MBhigh and MBlow elevation plant communities between study years. Mountain fens are typically smaller than boreal and subarctic peatlands due to strong topographic confinement in valleys and basins, and relatively small hydrological contributing areas (Patterson and Cooper 2007). Fens at lower elevations may be partially buffered from the adverse effects of climate change, depending on the size and geological nature of the watershed that support their hydrological regime. For example, fens in watersheds large enough to provide adequate groundwater flows to maintain shallow water tables during the growing season may not experience declines in NEP. It is important to note, however, that fen ecosystems within mountain watersheds, regardless of watershed contributing area, may experience increased water demands and drought stress due to climate change (Rood and others 2008) and human uses, such as ground water withdrawals (Cooper and others 2015). In such cases, the potential benefits of large watersheds may be diminished. The low rates of NEP at the low elevation sites were associated with warmer air temperatures, lesser snowpack, earlier melt-out, and lower mean growing season water tables than the high elevation sites. Coupled with predictions of reduced peak SWE and earlier snow melt associated with a warming climate, the results of this study suggest that the elevation ranges that provide climate conditions supportive of mountain fens in the Rocky Mountains may be narrowing, with their lower limits shifting upward in elevation. This supports our hypothesis that low elevation fens may be particularly susceptible to the effects of climate change. In cases where long-term groundwater and surface water base flows are unable to maintain shallow water tables in mountain fens, over time they may lose the organic soils that define them, through increased decomposition relative to production.

ACKNOWLEDGEMENTS We thank Drs. N.T. Hobbs and Phillip Chapman of Colorado State University for valuable assistance feedback on statistical analyses and model development. In addition, we thank the many students and field technicians whose contributions made this work possible. Funding for this study came from the United States Department of Agriculture––US Forest Service through the Rocky Mountain Research Station, in Fort Collins, Colorado. Data analyses and presentation were done using open-source R software including the rjags and ggplot2 package.

REFERENCES Adkinson AC, Syed KH, Flanagan LB. 2011. Contrasting responses of growing season ecosystem CO2 exchange to variation in temperature and water table depth in two peatlands in northern Alberta, Canada. J Geophys Res 116:1–17. Ashfaq M, Ghosh S, Kao SC, Bowling LC, Mote P, Touma D, Rauscher SA, Diffenbaugh NS. 2013. Near-term acceleration of hydroclimatic change in the western US. J Geophys Res Atmos 118:10,676–93. Aurela M. 2002. Annual CO2 balance of a subarctic fen in northern Europe: importance of the wintertime efflux. J Geophys Res 107:4607. Belyea LR. 2009. Nonlinear dynamics of peatlands and potential feedbacks on the climate system. In: Baird AJ, Belyea LR, Comas X, Reeve AS, Slater LD, Eds. Carbon Cycling in Northern Peatlands, Geophysical Monograph Series, vol. 184. Washington, D.C.: AGU. p 5–18. Blanken PD. 2014. The effect of winter drought on evaporation from a high-elevation wetland. J Geophys Res Biogeosciences 119:1354–69. Bubier JL, Bhatia G, Moore TR, Roulet NT, Lafleur PM. 2003. Spatial and temporal variability in growing-season net ecosystem carbon dioxide exchange at a large peatland in Ontario, Canada. Ecosystems 6:353–67. Cable JM, Ogle K, Barron-Gafford GA, Bentley LP, Cable WL, Scott RL, Williams DG, Huxman TE. 2013. Antecedent conditions influence soil respiration differences in shrub and grass patches. Ecosystems 16:1230–47. Chapin AFS, Woodwell GM, Randerson JT, Rastetter EB, Lovett GM, Baldocchi D, Clark DA, Harmon ME, Schimel DS, Valentini R, Wirth C, Aber JD, Cole J, Goulden ML, Harden JW, Heimann M, Howarth RW, Matson PA, Melillo JM, Mooney HA, Neff JC, Houghton RA, Pace ML, Ryan MG, Running W, Sala OE, Schlesinger WH, Schulze E, Chapin FS, Baldocchi DD, Harmon E, Cole JJ, Mcguire AD. 2006. Reconciling carbon-cycle concepts, terminology, and methods. Ecosystems 9:1041–50. Chimner RA, Lemly JM, Cooper DJ. 2010. Mountain fen distribution, types and restoration priorities, San Juan Mountains, Colorado, USA. Wetlands 30:763–71. Chimner RA, Cooper DJ, Parton WJ. 2002. Modeling carbon accumulation in Rocky Mountain fens. Wetlands 22:100–10. Chimner RA, Cooper DJ. 2003a. Influence of water table levels on CO2 emissions in a Colorado subalpine fen: an in situ microcosm study. Soil Biol Biochem 35:345–51.

D. J. Millar and others Chimner RA, Cooper DJ. 2003b. Carbon dynamics of pristine and hydrologically modified fens in the southern Rocky Mountains. Can J Bot 81:477–91. Christensen NS, Lettenmaier DP. 2007. A multimodel ensemble approach to assessment of climate change impacts on the hydrology and water resources of the Colorado River Basin. Hydrol Earth Syst Sci 11:1417–34. Christensen NS, Wood AW, Voisin N, Lettenmaier DP, Palmer RN. 2004. The effects of climate change on the hydrology and water resources of the Colorado River Basin. Clim Change 62:337–63. Cooper DJ, Andrus RE. 1994. Patterns of vegetation and water chemistry in peatlands of the west-central Wind River Range, Wyoming, USA. Can J Bot 72:1586–97. Cooper DJ. 1996. Water and soil chemistry, floristics, and phytosociology of the extreme rich High Creek fen, in South Park, Colorado, USA. Can J Bot 74:1801–11. Cooper DJ, Wolf EC, Ronayne MJ, Roche JW. 2015. Effects of groundwater pumping on the sustainability of a mountain wetland complex, Yosemite National Park, California. J Hydrol Reg Stud 3:87–105. Costigan KR, Bossert JE, Langley DL. 2000. Atmospheric/hydrologic models for the Rio Grande Basin: simulations of precipitation variability. Glob Planet Change 25:83–110. Davidson EA, Belk E, Boone RD. 1998. Soil water content and temperature as independent or confounded factors controlling soil respiration in a temperate mixed hardwood forest. Glob Change Biol 4:217–27. Davidson EA, Janssens IA, Luo Y. 2006. On the variability of respiration in terrestrial ecosystems: moving beyond Q10. Glob Change Biol 12:154–64. Frolking S, Roulet NT. 2007. Holocene radiative forcing impact of northern peatland carbon accumulation and methane emissions. Glob Change Biol 13:1079–88. Gelman A, Rubin DB. 1992. Inference from iterative simulation using multiple sequences. Statistical Sci 7:457–511. Heidel B, Jones G. 2006. Botanical and Ecological Characteristics of Fens in the Medicine Bow 676 Mountains, Medicine Bow National Forest Albany and Carbon Counties, Wyoming. Report 677 Prepared for: Medicine Bow-Routt National Forest. FS Agreement No. 02-CS-11020600-678 033 M8. Johnson CP, Pypker TG, Hribljan JA, Chimner RA. 2013. Open top chambers and infrared lamps: A comparison of heating efficacy and CO2/CH4 dynamics in a northern Michigan peatland. Ecosystems 16:736–48. Krajick K. 2004. All downhill from here ? Science 303:1600–2. Kruschke JK. 2013. Bayesian estimation supersedes the t test. J Exp Psychol Gen 142:573–603.

interactions, and forest carbon budgets. Proc Natl Acad Sci USA 108:9508–12. Millar, D. 2015. Climate controls on ecosystem-atmosphere carbon exchange and hydrological dynamics in Rocky Mountain fens. PhD dissertation. Graduate Degree Program in Ecology, Colorado State University, Fort Collins, CO. 107 pp. Moore C, Kampf S, Stone B, Richer E. 2014. A GIS-based method for defining snow zones: application to the western United States. Geocarto Int 30:62–81. Nijp JJ, Limpens J, Metselaar K, Peichl M, Nilsson MB, van der Zee SEATM, Berendse F. 2015. Rain events decrease boreal peatland net CO2 uptake through reduced light availability. Glob Change Biol 21:2309–20. Nilsson M, Sagerfors J, Buffam I, Laudon H, Eriksson T, Grelle A, Klemedtsson L, Weslien P, Lindroth A. 2008. Contemporary carbon accumulation in a boreal oligotrophic minerogenic mire––a significant sink after accounting for all C-fluxes. Glob Change Biol 14:2317–32. Parmesan C. 2006. Ecological and evolutionary responses to recent climate change. Annu Rev Ecol Evol Syst 37:637–69. Patterson L, Cooper DJ. 2007. The use of hydrological and ecological indicators for the restoration of drainage ditches and water diversions in a mountain fen, Cascade Range, California. Wetlands 27:290–304. ¨ quist M, Ottosson Lo¨fvenius M, Ilstedt U, Sagerfors J, Peichl M, O Grelle A, Lindroth A, Nilsson MB. 2014. A 12-year record reveals pre-growing season temperature and water table level threshold effects on the net carbon dioxide exchange in a boreal fen. Environ Res Lett 9:055006. Plummer M. 2011. Rjags: Bayesian graphical models using MCMC. http://CRAN.R-project.org/package=rjags. Ray AJ, Barsugli JJ, Averyt KB, Wolter K, Hoerling M, Doesken N, Udall B, Webb RS. 2008. Climate change in Colorado: A synthesis to support water resources management and adaptation. Colorado Water Conservation Board Rep., 52 pp. Regonda SK, Rajagopalan B, Clark M, Pitlick J. 2004. Seasonal cycle shifts in hydroclimatology over the western United States. J Clim 18:372–84. Ricker MC, Stolt MH, Zavada MS. 2014. Comparison of soil organic carbon dynamics in forested riparian wetlands and adjacent uplands. Soil Sci Soc Am J 78:1817. Riutta T, Laine J, Tuittila E-S. 2007. Sensitivity of CO2 exchange of fen ecosystem components to water level variation. Ecosystems 10:718–33. Rood SB, Pan J, Gill KM, Franks CG, Samuelson GM, Shepherd A. 2008. Declining summer flows of Rocky Mountain rivers: Changing seasonal hydrology and probable impacts on floodplain forests. J Hydrol 349:397–410.

Maanavilja L, Riutta T, Aurela M, Pulkkinen M, Laurila T, Tuittila E-S. 2011. Spatial variation in CO2 exchange at a northern aapa mire. Biogeochemistry 104:325–45.

Rydin H, Jeglum J. 2006. The biology of peatlands. Oxford University Press. p 239. Sagerfors J, Lindroth A, Grelle A, Klemedtsson L, Weslien P, Nilsson MB. 2008. Annual CO2 exchange between a nutrientpoor, minerotrophic, boreal mire and the atmosphere. J Geophys Res Biogeosciences 113:1–15.

Mast MA, Wickland KP, Striegl RT, Clow DW. 1998. Winter fluxes of CO2 and CH4 from subalpine soils in Rocky Mountain National Park, Colorado ecosystems. Glob Biogeochem Cycles 12:607–20.

Sanchez A, Hughes NM, Smith WK. 2015. Importance of natural cloud regimes to ecophysiology in the alpine species, Caltha leptosepala and Arnica parryi, Snowy Range Mountains, southeast Wyoming, USA. Funct Plant Biol 42:186–97.

Melillo JM, Butler S, Johnson J, Mohan J, Steudler P, Lux H, Burrows E, Bowles F, Smith R, Scott L, Vario C, Hill T, Burton A, Zhou YM, Tang J. 2011. Soil warming, carbon-nitrogen

Schimelpfenig DW, Cooper DJ, Chimner RA. 2014. Effectiveness of ditch blockage for restoring hydrologic and soil processes in mountain peatlands. Restor Ecol. 2:257–65.

Letts MG, Roulet NT, Comer NT, Skarupa MR, Verseghy DL. 2000. Parametrization of peatland hydraulic properties for the Canadian land surface scheme. Atmos Ocean 38:141–60.

Mountain Peatlands Range from CO2 Sinks at High Elevations to Sources at Low Elevations Schneider J, Kutzbach L, Wilmking M. 2012. Carbon dioxide exchange fluxes of a boreal peatland over a complete growing season, Komi Republic, NW Russia. Biogeochemistry 111:485–513. Sullivan PF, Arens SJT, Chimner RA, Welker JM. 2007. Temperature and microtopography interact to control carbon cycling in a high arctic fen. Ecosystems 11:61–76. Tesoriero AJ, Terziotti S, Abrams DB. 2015. Predicting redox conditions in groundwater at a regional scale. Environ Sci Technol 49:9657–64. Tufford DL. 2011. Shallow water table response to seasonal and interannual climate variability. Trans ASABE 54:2079–86. Van Meter KJ, Basu NB. 2015. Catchment legacies and time lags: a parsimonious watershed model to predict the effects of legacy storage on nitrogen export. PLoS One 10: e0125971. Vidon P. 2012. Towards a better understanding of riparian zone water table response to precipitation: surface water infiltration, hillslope contribution or pressure wave processes? Hydrol Process 26:3207–15. Walker TN, Garnett MH, Ward SE, Oakley S, Bardgett RD, Ostle NJ. 2016. Vascular plants promote ancient peatland carbon loss with climate warming. Glob Change Biol:n/a–n/a.

Warshall P. 1994. The Madrean Sky Island Archipelago: A Planetary Overview. In: Biodiversity and Management of the Madrean Archipelago: The Sky Islands of Southwestern United States and Northwestern Mexico (ed De Bano LF, Folliott PF O-RA), pp. 6–18. US Department of Agriculture, Fort Collins, Colorado. Wickland KP, Striegl G, Mast MA, Clow DW. 2001. Carbon gas exchange at a southern Rocky Mountain wetland, 1996–1998. Glob Biogeochem Cycles 15:321–35. Wohlfahrt G, Cremonese E, Hammerle A, Ho¨rtnagl L, Galvagno M, Gianelle D, Marcolla B, di Cella UM. 2013. Tradeoffs between global warming and day length on the start of the carbon uptake period in seasonally cold ecosystems. Geophys Res Lett 40:6136–42. Wu J, Roulet NT, Sagerfors J, Nilsson MB. 2013. Simulation of six years of carbon fluxes for a sedge-dominated oligotrophic minerogenic peatland in Northern Sweden using the McGill Wetland Model (MWM). J Geophys Res Biogeosciences 118:795–807. Yurova A, Wolf A, Sagerfors J, Nilsson M. 2007. Variations in net ecosystem exchange of carbon dioxide in a boreal mire: Modeling mechanisms linked to water table position. J Geophys Res Biogeosciences 112:1–13.