This item is only available as the following downloads:
1 EFFECTS OF PERMAFROST THAW ON TUNDRA CARBON BALANCE OVER SPACE AND TIME By ELIZABETH FAY BELSHE A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2013
2 2013 Elizabeth Fay Belshe
3 To my friends Pep Sanjaun, Allie Shenkin, Peter Kelly, Paulo Brando, Caitlin Hicks Pries, Christina Feilding Jessica Diller
4 ACKNOWLEDGMENTS I would like to thank Dr. Ben Bolker for the countless hours of mentoring and support he provided on ecological modeling, which ultimately gave me the confidence and knowledge to understand and analyze the data in my dissertation. Next, I would like to thank Dr. Ted Schuur for his guidance, advice, and endless support throughout my dissertation. I also want to thank Dr. Rosvel Bracho for his help with site level data analysis and helpful comments and guidance and thanks to Caitlin Hicks Pries and Pep Sanjuan for allowing me to think out loud and work through my dissertation Thanks to Forrest Stevens for support with spatial analysis and ArcGIS, Alexander Shenkin for computer and moral support, and Paulo Brand o for R support. For valuable field support and ideas during the initial set up of eddy covariance tower I thank Christian Trucco and Dr. Jason Vogel. Thanks to Sasha Ivans at Campbell Scientific for eddy covariance technical support and training. I would also like to thank UNAVCO for providing GPS equipment, training, and support, which is supported by the National Science Foundation and National Aeronautics and Space Administration under NSF Cooperative Agreement No. EAG 0735156. Finally, I would like to thank the Department of Energy, Graduate Research Environmental Fellowship for support during the last three years of my dissertation.
5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ ............... 4 LIST OF TABLES ................................ ................................ ................................ ........................... 7 LIST OF FIGURES ................................ ................................ ................................ ......................... 8 ABSTRACT ................................ ................................ ................................ ................................ ..... 9 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .................. 11 Permaf rost ................................ ................................ ................................ ............................... 11 Permafrost and Carbon ................................ ................................ ................................ ........... 12 Future of Permafrost ................................ ................................ ................................ ............... 13 Changes in Permafrost and Implications for the Global Carbon Cycle ................................ .. 14 2 QUANTIFICATION OF UPLAND THERMOKARST FEATURES WITH HIGH RESOLUTION REMOTE SENSING ................................ ................................ .................... 15 Background ................................ ................................ ................................ ............................. 15 Material and Methods ................................ ................................ ................................ ............. 17 Site Description ................................ ................................ ................................ ............... 17 Land Cover Classification ................................ ................................ ............................... 18 Im age acquisition and field sampling ................................ ................................ ....... 18 Data preparation ................................ ................................ ................................ ....... 20 Supervised classification ................................ ................................ .......................... 21 Cl assification accuracy assessment ................................ ................................ .......... 21 Results ................................ ................................ ................................ ................................ ..... 22 Discussion ................................ ................................ ................................ ............................... 24 Concluding Remarks ................................ ................................ ................................ .............. 27 3 INCORPORATING SPATIAL HETEROGENEITY CREATED BY PERMAFROST THAW INTO A LANDSCAPE CARBON ESTIMATE ................................ ....................... 36 Background ................................ ................................ ................................ ............................. 36 Material and M ethods ................................ ................................ ................................ ............. 38 Site Description ................................ ................................ ................................ ............... 38 Eddy Covariance Measurements ................................ ................................ ..................... 39 Eddy Covariance Data Handling ................................ ................................ ..................... 40 Environmental Measurements ................................ ................................ ......................... 41 Landscape Properties ................................ ................................ ................................ ....... 41 Estimation of Landscape Scale Carbon Exchange ................................ .......................... 44 Gap filling strategy 1: generalized additive models ................................ ................. 45
6 Gap filling strategy 2: non linear regressio ns ................................ .......................... 47 Model Performance ................................ ................................ ................................ ......... 48 Results ................................ ................................ ................................ ................................ ..... 49 Landscape Heterogeneity of the EC Footprint ................................ ................................ 49 Model Performance ................................ ................................ ................................ ......... 50 Predictions of Ecosystem C Balance ................................ ................................ ............... 50 Predictions of Landscape Heterogeneity of C Flux ................................ ......................... 51 Discussion ................................ ................................ ................................ ............................... 52 Quantifying Spatial Heterogeneity ................................ ................................ .................. 52 Incorporating Spatial Heterogeneity into EC C Estimates ................................ .............. 54 Effects of Spatial Heterogeneity on C Flux ................................ ................................ ..... 55 Concluding Remarks ................................ ................................ ................................ .............. 57 4 TUNDRA ECOSYSTEMS OBSERVED TO BE CO2 SOURCES DUE TO DIFFERENTIAL AMPLIFICATION OF THE CARBON CYCLE ................................ ..... 66 Background ................................ ................................ ................................ ............................. 66 Material and Methods ................................ ................................ ................................ ............. 68 Literature Search ................................ ................................ ................................ ............. 68 Data Analysis ................................ ................................ ................................ ................... 70 Results ................................ ................................ ................................ ................................ ..... 72 Changes in CO 2 Flux Over Time ................................ ................................ ..................... 72 Relationships Between CO 2 Flux and Spatial Differences in Climate ............................ 74 Discussion ................................ ................................ ................................ ............................... 75 5 CONCLUSION ................................ ................................ ................................ ....................... 85 LIST OF REFERENCES ................................ ................................ ................................ ............... 89 BIOGRAPHICAL SKETCH ................................ ................................ ................................ ....... 105
7 LIST OF TABLES Table page 2-1 Characterization of land cover classes. .............................................................................. 29 2-2 Spatial extent and percent of total area covered by each class within the study area. .......30 2-3 Summary of supervised classification accuracy assessment. ............................................31 3-1 Model preformance of GAM and NL models....................................................................60 3-2 Carbon estimates for the wind directions. ..........................................................................60
8 LIST OF FIGURES Figure page 2 1 IKONOS image of the EML watershed and adjacent hill slope. ................................ ....... 32 2 2 Mapped supervised classification of the EML watershed. ................................ ................ 33 2 3 So il organic layer thickness of land cover classes ................................ ............................. 34 2 4 Max imum seasonal thaw depth of land cover classes ................................ ...................... 35 3 1 Map of microtopography and 360 tran sects surrounding the EC tower. ........................... 61 3 2 R elationship s betw een microtopography and ecosystem properties ................................ 62 3 3 Relationships between r oughness and net ecosystem exchange. ................................ ....... 63 3 4 Carbon fluxes for the three methods of prediction. ................................ ........................... 64 3 5 Carbon fluxes for different wind direction s ................................ ................................ ...... 65 4 1 Map of the Wo rld showing the 32 sites included in this meta analysis. ............................ 80 4 2 Temporal trends of net ecosystem CO 2 exchange. ................................ ............................ 81 4 3 Estimates of temporal trends of net ecosystem exchange. ................................ ................. 82 4 4 Linear relationships between mean annual temperature and carbon fluxes. ..................... 83 4 5 Estimates of re lationships between mean annual temperature carbon fluxes. ................... 84
9 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment o f the Requirements for the Degree of Doctor of Philosophy EFFECTS OF PERMAFROST THAW ON TUNDRA CARBON BALANCE OVER SPACE AND TIME By Elizabeth Fay Belshe May 2013 Chair: Edward Schuur Major: Botany Climate induced changes to permafrost are altering high latitude landscapes in ways that could increase the vulnerability of the vast soil carbon pools of the region. T hawing of permafrost is a heterogeneous process t hat creates a complex landscape. In addition to the thickening of soil active layers, localized thermokarst features form when ice rich permafrost thaws and the ground subsides. By p er forming a supervised classification on a high resolution IKONOS image, I detected and mapped small, irregular thermokarst features occu rring within an upland watershed in Interior Alaska Twelve percent of the Eight Mile Lake (EML) watershed has undergone thaw and thermokarst predominantly form in valleys where tussock tundra resides. By incorporating spatial variation created by permafro st thaw into an eddy covariance carbon estimate, I estimated the carbon balance o f tundra containing thermokarst During the two 6 month periods measured permafrost thaw increased the amplitude of the carbon cycle by stimulating both carbon release and sequestration, while the ecosystem remained a carbon sink at the landscape scale To put these results into context and determine if tundra ecosys tem s were carbon sources or sinks I compile d 40 years of CO 2 flux observations from 54 studies spanning 32 sites across northern high latitudes Using time series analysis, I investigated if seasonal or annual CO 2
10 fluxes have changed over time, and whether spatial differences in mean annual temperature could help explain temporal changes in CO 2 flux Growing season net CO 2 uptake has definitely increased since the 1990s; the data also suggest (albeit less definitively) an increase in winter CO 2 emissions, especially in the last decade. In spite of the uncertainty in the winter trend, we estimate that tundra sites were annual CO 2 sources from the mid 1980s until the 2000s; and data from the last seven years show that tundra continue s to emit CO 2 annually. CO 2 emissions exceed CO 2 uptake across the range of temperatures that occur in the tundra biome. Taken together, these data suggest that despite increases in growing season uptake, tundra ecosystems are currently CO 2 sources on an annual basis.
11 CHAPTER 1 INTRODUCTION Permafrost A defining characteristic of many high latitude ecosystems is the presence of permafrost, whi ch occurs within 24 % of the land area in the northern hemisphere ( Anisimov & Nelson 1996 ; Brown et al. 1998 ) Permafrost is defined as any earth mate rial that remains at or below 0C for two or more consecutive years. Since p ermafrost ground is determined by tempe rature and not substrate material, i t can encompass a wide range of material from solid rock, to mineral and organic soil, to ice ( Davis 2001 ) The surface zone of permafrost that thaws during the summer and completely refreezes during the winter is defined as the active layer ( Shur et al. 2005 ) This active layer ranges from a few tens of centimeters to more than several meters thick ( Davis 2001 ) Although active layer thickness is largely controlled by regional climate, other factors such as topographical aspect, vegetation su bstrate, and organic matter insulation, and heat conductance of the overlying soil also play a role in determining the extent of seasonal thaw. Global permafrost distribution falls into two broad zones dictated by climate: the zone of continuous permafrost where climate is favorable to permafrost formation and stability and the zone of discontinuous permafrost where climate is neutral to permafrost and it can be either present or absent depending on local conditions ( Shur & Jorgenson 2007 ) Because climate is dynamic, the distribution of permafrost is continuously changing to be in equilibrium with contemporary climate. While climate is the dominate factor th at determines the geographic extent of permafrost, ecosystem level processes such as soil development, vegetation succession, and disturbance interact to affect the stability and distribution of permafrost ( Shur & Jorgenson 2007 ) Permafrost that forms under favorable conditions consists of a continuous subsurface layer overlain by a
12 very thick active layer. Ecosystem level processes, such as vegetation changes and soil development, can then act upon and modify the upper permafrost. During ecosystem succession, vegetation cover increases and organic matter accumulates, which in turn buffers the permafrost from air temperatures and results in the upward movement of the permafrost table. This thinning of the active layer is accompan ied by an accumulation of ice, which forms an ice rich intermediate layer between th e active layer and permafrost. This ecosystem modified permafrost is more thermally stable because of the protective organic layer but more susceptible to thaw because of its high ice content ( Shur & Jorgenson 2007 ) Permafrost and Carbon Carbon is incorporated into permafrost by mechanisms that redistribute it away from the surface to d eeper in the soil profile. Organic carbon is initially added to the soil by primary producers and factors such as rooting de pth of plants determine the initial depth distribution of new carbon inputs ( Schuur et al. 2008 ) Then i nputs to the soil from either organic matter deposition or various sedimentation processes cause soil surface accretion. This causes a vertical increase in the soil surface and an upward movement of the permafrost table which results in organic matter at the bottom of the active layer becoming incorporated into the perm afrost ( Schuur et al. 2008 ) In addition cryoturbation (physical mixing of soil horizons that occurs during freeze thaw processes) helps to move organic matter down into the permafrost ( Michaelson et al. 1996 ) These processes result in a variety of organic matter s ubstrates (from labile to recalcitrant) becoming incorporated into permafrost. Then many factors including temperature, hydrology, and substrate quality interact to stabilize carbon in these soils ( Hobbie et al. 2000 ) T he total permafrost soil carbon pool in the northern circumpolar zone is estimated to be 1672 Pg, which is twice the size of the atmospheric carbon pool ( Schuur et al. 2008 ; Tarnocai et al. 2009 )
13 Future of Permafrost Global s urface air temperature has increased 0.35 C decade 1 ( Polyakov et al. 2002 ; Hinzman et al. 2005 ; Serreze & Francis 2006 ; Euskirchen et al. 200 7 ) and temperatures are predicted to rise by up to 6.5 C in northern high latitudes by the year 2100 ( IPCC 2007 ) It is estimated that 25 90% of areas with near surface permafrost will transition to seasonally frozen ground by the year 2100 ( Anisimov & Nel son 1996 ; Saito et al. 2007 ; Lawrence et al. 2008 ; Jafarov et al. 2012 ) In add ition, the rate and areal extent of permafrost thaw has already begun to increase ( Camill 2005 ; Hinzman et al. 2005 ; Jorgenson et al. 2006 ; Schuur et al. 2008 ; Osterkamp et al. 2009 ; Romanovsky et al. 2010 ) and active layer thickening has been observed ( Zhang et al. 2005 ) As mean near surf geographic extent of permafrost and active layer thickness are expected to change ( Anisi mov et al. 2007 ) But i t is difficult to infer changes in permafrost from changes in air temperature b ecause various disturbances and feedbacks with hydrology, vegetation, and geomorphology complicate the response of permafrost to warming ( Stieglitz et al. 2003 ; Grosse et al. 2011 ) And in addition to regional level changes in permafrost distribution and active layer thickness, there can also be localize d permafrost degradation that results in th ermokarst formation. Thermokarst topography forms when ice rich permafrost thaws, ground ice melts, surface and subsurface drainage occurs and the ground surface subsides due to volume loss during ice water phase transition ( Hinzman et al. 2005 ) The amount of settlement is dependent on the amount and types of ice that occur in the permafrost, which in turn is dependent on the slope, soil texture, hydrolog y, vegetation and climatic and geologic history of the area. So, a t the landscape scale there can be a wide diversity of thermokarst features because of differences in local conditions ( Jorgenson & Osterkamp 2005 ) Because of these complex interactions and the diversity of
14 permafrost composition the predictability of permafrost thaw and thermokarst formation i s a formidable challenge Changes in Permafrost and Implications for the Global Carbon Cycle Permafrost thaw could result in a strong positive feedback to climate change if a portion of the estimated 1670 Pg of soil organic carbon stored in permafrost regions is released to the atmosphere ( Schuur et a l. 2008 ; Tarnocai et al. 2009 ) The future carbon balance of high latitude ecosystems is dependent on the sensitivity of biological processes (photosynthesis and respiration) to the physical changes occurring with permafrost thawing. Higher soil temperatures increase active layer thickness and create the rmokarst terrain, which can change soil hydrology and in turn can either increase or decreas e gross primary production and ecosystem respiration. T he balance between these large opposing fluxes will depend on the cumulative effect of the landscape level ph ysical changes that occur with permafrost thaw Therefore, to understand the effects of climate change on the rate and magnitude of carbon exchange, we must quantify both the effects of permafrost thaw on carbon fluxes and the spatial distribution of thaw features, such as thermokarst, across the landscape. The objectives of this dissertation are to: 1) develop methodology to detect and classify typical thermokarst features occurring in an upland tundra watershed, 2) quantify the effect of permafrost thaw and thermokarst on carbon exchange of an upland tundra ecosystem ; and 3) using meta analysis assess whether seasonal or annual carbon fluxes within the tundra biome have changed over the past four decades as a result of changes already occurring in high latitude regions
15 CHAPTER 2 QUANTIFICATION OF UPLAND THERMOKARST FEATURES WITH HIGH RESOLUTION REMOTE SENSING Background Temperature is increasing in northern high latitudes in response to the radiative forcing associated with increasing gre enhouse gases and changes in albedo ( Chapin et al. 2005 ) Surface air temperature has increased 0.35 C decade 1 ( Polyakov et al. 2002 ; Hinzman et al. 2005 ; Serreze & Francis 2006 ; Euskirchen et al. 2007 ) and temperatures are predicted to rise by up to 6.5 C in northern high latitudes by the year 2100 ( IPCC 2007 ) This warming affects the spatial distribution of permafrost, which currently occurs within 24% of the ice free land area in the northern hemisphere ( Zhang et al. 1999 ) It is estimated that 25 90% of areas with near surface permafrost will transition to seasonally frozen ground by the year 2100 ( Anisimov & Nelson 1996 ; Saito et al. 2007 ; Lawrence et al. 2008 ; Jafarov et al. 2012 ) and rate and areal extent of permafrost thaw has already begun to increase ( C amill 2005 ; Hinzman et al. 2005 ; Jorgenson et al. 2006 ; Schuur et al. 2008 ; Osterkamp et al. 2009 ; Romanovsky et al. 2010 ) This thaw could result in a strong positive feedback to climate change if a portion of the estimated 1670 Pg of soil organic carbon in permafrost regions is released to the atmosphere ( Schuur et al. 2008 ; Tarnocai et al. 2009 ) Understanding the relationships among warming, permafrost thaw, and carbon cycling is a challenge because permafrost thaw is temporally dynamic and spatially heterogeneous. Various disturbances and feedbacks with hydrology, vegetation, and geomorphology additionally complicate the response of permafrost to warming ( Grosse et al. 2011 ) In addition to the Belshe EF, EAG Schuur G Grosse (submitted) Quantification of upland thermokarst features with high resolution remote sensing. Environmental Research Letters
16 predicted thickening of the active layer (seasonally thawed surface layer) ( Anisimov et al. 2007 ) thawing of ice rich permafrost can result in thermokarst formation ( Hinzman et al. 2005 ; Jorgenson & Osterkamp 2005 ) Thermokarst landforms often have very different physical and biogeochemical environments compared to the surrounding un affected permafrost landscape, with impacts on both rates of carbon uptake and emission ( Walter et al. 2007 ; Vogel et al. 2009 ; Lee et al. 2010 ; Belshe et al. 2012 ; Jones et al. 2012 ) Ultimately, net landscape scale carbon exchange will depend on the cumulative response of the mosaicked landscape created by permafrost thaw. Therefore, to estimate the rate and magnitude of carbon efflux from the landscape it is imperative to quantify the size and distribution of thermokarst landforms. The wide diversity of thermokarst features depends on local permafrost and ground ice conditions, landscape position, hydrology, soil texture, and surrounding vegetation ( Jorgenson & Osterkamp 2005 ) Sixteen primary modes of surface permafrost degradation have been identified, ranging from features with substantial thaw settlement (2 20 m) such as glacial thermokarst and drained lake basins to features with relatively low subsidence ( 0.2 2 m) such as thermokarst gullies and non patterned thawed ground (Jorgenson & Osterkamp 2005). One of the major limitations for estimating the size, number, and spatial distribution of thermokarst features is attaining adequate spatial coverage because high latitude ecosystems are vast, isolated, and often inaccessible. Remote se nsing allows us to observe these areas in a cost effective manner and on the timescales needed to answer questions about the rapidly changing landscape. So far, the remote detection of thermokarst has been limited to relatively large and visually distinct features, such as thermokarst lakes ( Jones et al. 2011 ; Morgenstern et al. 2011 ) drained lake basins ( Hinkel et al. 2003 ; Wang et al. 2012 ) retrogressive thaw slumps ( Lantuit & Pollard 2008 ; La ntuit et al. 2012 ) and thermo erosion gullies ( Godin & Fortier 2012 ) However,
17 to capture the true magnitude of thermokarst, it is also necessary to expand o bservation capabili ties towards detection of small, subtle thermokarst forms. With recent advancements in spectral and spatial resolution this is becoming possible. The objective of this study is to assess the feasibility of very high resolution remote sensing for detecting small irregular thermokarst features (pits, gullies, and long water tracks) in an upland permafrost area in the discontinuous permafrost zone where such features are numerous. Similar thermokarst gullies and pits occur throughout the ~64,000 ha Toklat Basi n that extends 80 km to the west of our study site in Interior Alaska, and in many other high latitude upland areas. By performing a land cover classification targeting thermokarst features in IKONOS satellite imagery, we quantified the spatial extent of t hermokarst and determine their spatial distribution and landscape position. Using this information, along with field observations of thaw depth and organic soil depth, we exemplarily assess the vulnerability of soil carbon pools in a typical upland area in the discontinuous permafrost zone of Interior Alaska. Material and Methods Site Description adjacent hill slopes located in the northern foothills of the Alaska Range n ear Denali National Park and Preserve, Interior Alaska ( Schuur et al. 2007 ; Schuur et al. 2009 ) The study areas covers 10.4 km 2 and is located at ~700m elevation on a gently sloping (~5%), north facing glacial terminal moraine that dates back to the Early Pleistocene. Nenana gravel constitutes the major bedrock formation of the area and due to multiple glacial advances and retreats a landscape of rolling hills and valleys has formed ( Wahrhaftig 1958 ) The hilltops have gravelly soils associated with glacial til l and are covered by a thin aeolian silt cap, while the lower slopes
18 and valleys have thick peat accumulations over fine grained soils associated with colluvial deposits ( Osterkamp et al. 2009 ) The study area is underlain by permafrost but it is a lso mapped within the vulnerable band of discontinuous permafrost near the point of thaw due to the combination of its elevation and geographic position ( Yocum et al. 2006 ; Romanovsky et al. 2007 ) Permafrost at EML is classified as climate driven, ecosystem protected ( Shur & Jorgenson 2007 ) because the permafrost was formed under an earlier cold climate independent of ecological conditions and later protected from degradation by vegetation and soil d evelopment ( Osterkamp et al. 2009 ) The long term mean annual air temperature (1976 2009) of the area is 1.0 C and the growing season mean temperature is 11.2 C, with monthly averages ranging from 16 C in D ecember to +15 C in July. The long te rm annual precipitation is 378 mm with a growing season precipitation of 245 mm (National Climatic Data Center, NOAA). Deep permafrost temperature has been measured in a borehole at the site since 1985 and during this time mean annual ground temperature ro se by 1C (1989 1998) and thermokarst terrain has developed and expanded ( Osterk amp et al. 2009 ) Land Cover Classification Image acquisition and field sampling A multi spectral IKONOS satellite image covering the EML watershed and adjacent hillslope was collected on the 25 th June 2008 (Figure 2 1). IKONOS Geo 1 products have a spatial resolution of 3.2 meters for multispectral bands with a position accuracy of 53 meters. To improve this georectification, eleven ground control points collecte d in the field with a differential global positioning system (DGPS) along with nearest neighbor resampling and a 2 nd order polynomial transformation were used to rectify the image, achieving an average root mean square error ( RMSE ) of 0.98 meters.
19 Informa tion on vegetation and soil was collected within 8 a priori selected land cover types ( 2 1). An unsupervised classification of the IKONOS image was used to aid in selection of these potential classes along with field knowledge of the study area. Our sampli ng strategy focused on spatial coverage, therefore rapid assessment quantification metrics were used. Approximately 10 field training sites were sampled within each class, resulting in a total of 78 training sites distributed through out the watershed (Tabl e 2 1, Figure 2 1). At each site, four contiguous (2 x2 m) plots were established within homogeneous areas that were representative of was estimated and geog raphic location was taken with DGPS. Vegetation was further categorized by the percent cover of functional groups (graminoid, herbaceous, non vascular, shrub, tree) and percent cover of each species. Due to identification difficulties and time limitations, Salix spp., Alnus spp., Dicranum spp. and Sp h agnum spp. were only identified to genera and rare non vascular plants were grouped into a non These data were used to further characterize the eight land cover types (table 1), which we compared to the Denali landcover mapping project (DLMP; Stephens et al, 2001) and the Circumpolar Arctic vegetation map (CAVM; Walker et al, 2005). At the end of the growing season (late August) active layer thickness (depth to permafrost) was measured at each site by inserting a metal rod until the impermeable permafrost layer was reached. Differences in thaw depth among land cover classes were evaluated using an ANOVA with a Tukey post hoc test and analysis was p er for med in R (R Core Development Team 2011). In addition, depth of the near surface organic horizon was measured in small soil pits dug at each site. Because deep organic horizons (40 60cm) occurred in some classes and due to field time limitations, pits were only dug to ~20 cm.
20 If the organic horizon extended below 20 cm the site was considered highly organic and grouped into a >20 cm category. Fine scale differences in elevation of the study area were measured with a survey grade DGPS. One GPS unit (Trimble 5400) was placed at a nearby USGS geodetic marker (WGS84, unit (Trimble 5400) secured to a backpack, a kinematic survey was conducted by walking transects of the 10.4 km study site. Geographic position and elevation were collected at 5 s intervals, yielding a total of 68,000 data points. These data were post processed with methodology developed by UNAVCO using Trimble Geomatics Office software (Dayton, Ohio). Data prep aration A combination of DGPS point measurements of elevation and USGS Hydrologic Line Graphs Maps were used to create a digital elevation model (DEM) with 5 m grid cell size of the study area with the Topo to Raster function in the 3D Analyst Tool of ArcG IS 9.3. Slope was derived from the resulting DEM after it was smoothed using mean focal statistics (10 x 10 cell rectangle) with the Spatial Analyst tool in GIS. IKONOS bands 3 and 4 were used to calculate the normalized difference vegetation index (NDVI), which is a measure of photosynthetically ( Raynolds et al. 2006 ) We expected that NDVI could potentially help separate thermokarst classes because vegetation within thermokarst at our study site has been shown to have higher productivity than surrounding tussock tundra ( Schuur et al. 2007 ) To reduce the spectral overlap between some training classes, we also performed a principal components analysis (PCA) on IKONOS bands 1 4. Based on visual inspection of spectral separation, we used the first three PCA bands in subsequent analysis.
21 Supervised classification For the final classificat ion we combined six data layers (PCA1, PCA2, PCA3, NDVI, elevation, slope) and performed a maximum likelihood classification based on spectral signatures of the training sites. The spectral separability of individual training data were visually evaluated i n scatter plots combining the different data layers. Classes were merged if their spectra appeared to have insufficient separation. Originally the spruce, alder and willow classes were evaluated independently, but because they spectrally were very similar they were merged into a mixed riparian class. Post classification majority filtering was done to remove small numbers of pixels within larger contiguous classes and to achieve a smoother final classification. All classification analysis was done with ERDAS Imagine 9.1 software (Atlanta, GA). Classification accuracy assessment Independent validation sites representative of each land cover class were collected during the field campaign in 2009 (78 sites total; ~10 within each class) and from high resolution true color aerial images taken of the watershed in 2003 (203 total). Cros s tabulated error matrices of the mapped classes and validation sites were used to assess classification accuracy ( Congalton & Green 1999 ) ates the probability of a reference pixel being correctly accurately represents that category on the ground (commission error). The Kappa statistic incorporat es the off diagonal elements of the error matrices (i.e. classification errors), which represents agreement obtained after removing the proportion that could be expected by chance ( Jensen 2005 )
22 Results The EML watershed grades in elevation from 680 m in the south to 632 m in the north, where the lake is situated. Superimposed upon this overall gradient in elevation are rolling hills and valleys that differ in elevation by roughly 15 m. The final thematic map derived from the supervised classification shows that hil led areas correspond with the Hilltop Shrub (HS) class, which transition to the tussock tundra (T T) class found in the valleys (F igure 2 2). The two classes associated with permafrost degradation ( Thermokarst, TK and Thermokarst Pool, TP) mainly occur in the valleys of the watershed within the tussock tundra class ( F igure 2 2). In general, thermokarst features at EML are e long ated (~ 300 700 m) slightly incised (<1 m) gullies that drain from the valleys towards the lake (see inset of F igure 2 2 for a clo se up of typical thermokarst patterns). On average, ground subsidence within thermokarst features is less than 1 meter, and while individual gullies are relatively narrow (1 5 m wide), aggregated together they can affect areas up to 75 m wide (F igure 2 2). The mapped thermokarst classes combined (TK and TP) cover 1.31 km 2 which is 12.1 % of the 10.4 km 2 study area (T able 2 2). Tussock tundra (TT) covers 2.4 km 2 (22.9% of the total area), Hilltop shrub (HS) covers 1.7 km 2 (15.8%), Shrub Tundra Transition (S T) covers 2.3 km 2 (22%), and the Mixed Riparian (MP) covers 2.2 km 2 (21%; T able 2 2). Assuming that the thermokarst areas originally were all covered with tussock tundra, we estimate that the unthawed tussock tundra originally covered 3.7 km2 and by 2008 approximately 1/3 (34.6%) had undergone permafrost thaw. All training sites of thermokarst classes (TK, TP) and tus sock tundra (TT) had soil organic layers greater than 20 cm thickness, while all other classes had organic layers shallower than 20 c m (Fig ure 2 3). Our field method of categorizing all training sites with organic layers >20 cm as highly organic does not a llow us to statistically differentiate between classes because we did not measure variance within sites with thick organic layers due to the challenge of
23 sampling into frozen soil with limited field time However, previous work in the same watershed report ed that soil organic layers range from 373.1 cm at minimally thawed tundra sites to 544.8 cm at extensively thawed sites ( Hicks Pries et al. 2012 ) These site specific data along with our field data across the watershed show that permafrost subsidence at EML is occurring in areas with relatively large soil organic layers. Thermokarst classes (TK, TP) have greater maxim um seasonal thaw depths than other land cover classes (F igure 2 4). On average, thaw depth within TK and TP classes was 73 and 85 cm, respectively. The Mixed Riparian class (RM) did have a few sites near the lake inlet that exhibited deep thaw depths (66 & 81 cm), but on average all other land cover classes had thaw depths of 44 cm or lower. RM sites near the lake inlet most likely exhibited deeper thaw depths due to the thermal influence of a thaw bulb developing beneath the lake. Notably, the tussock tund ra (TT) class had much lower and less variable thaw depths than the thermokarst classes, even though they are adjacent to each other. This confirms that soils within thermokarst are seasonally thawing to a greater extent than any other land cover class. Th e overall classification accuracy of the IKONOS image was 83%, and there was igh, ranging from 6 9% to 100% (T able 2 accuracy of 75% and 100%, for TK and TP classes respectively. For the thermokarst classes, the highest confusion occurred between the TK and ST class, with 15 % of TK training pixels being classified as ST. This misclassification most likely occurred because sloped sites within the TK class have very similar vegetation compositions as the ST class. Both are dominated by short stature shrubs and have similar spec ies compositions (T able 2 1). Previous work at EML has
24 shown that the vegetation composition shifts from tussock forming sedges to shrubs as thermokarst develop ( Schuur et al. 2007 ) This is similar to the vegetation change that occurs at the transition between tussock tundra valleys and shrubby hilltops. The other obvi ous misclassification occurred between the MR and HS classes, with 31% of MR tra ining pixels classified as HS (T able 2 3). We believe this confusion occurred because the MR class is a mosaic of vegetation types that include tall shrubs that are spectrally similar to the medium stature shrubs occurring on the hilltops. Additionally, the classification was not well resolved in the southern portion of the watershed and overestimates the homogeneity and dominance of the RM class at the southern tip. This incong ruity was caused by the low number of training sites in this steeply sloped region. However we achieved an overall acceptable accuracy, with thermokarst classes having comparable accuracy to other classes in the landscape. Discussion Twelve percent of the landscape within EML watershed has undergone permafrost thaw resul ting in thermokarst formation (T able 2 2). Although thermokarst features are widespread throughout the watershed, they predominately occur within the low lying valleys dominated by tussock tundra (F igure 2 2). This pattern of ground subsidence is dictated by the initial amount and distribution of ice rich permafrost, which in turn resulted from the geology and soils of the area. The EML watershed experienced multiple glacial advances and ret reats that created a landscape of rolling hills and valleys ( Wahrhaftig 1958 ) The hills have gravelly soils consisting of glacial till covered by a thin aeolian silt cap with low amounts of ice within the permaf rost. Alternatively, the lower slopes and valleys have fine grained soils associated with colluvial deposits and ice rich horizon s containing up to 50% visible ice content ( Osterkamp et al. 2009 ) These ice layers and lenses formed over time as pea t accumulated in the valleys and water was trapped by upward freezing of the permafrost surface ( Osterkamp et al. 2009 ) Although the peat
25 accumulations initially acted to build and protect the permafrost, the co ncurrent increase in ice ultimately made permafrost in the valleys vulnerable to thermokarst formation. Through continued warming, insulation of soil by snow accumulating in subsided areas, and thermal erosion (from water being redistributed from higher to lower microtopographical areas) permafrost thaw has increased and thermokarst gullies and pools have expanded ( Kane et al. 2001 ; Osterkamp et al. 2009 ) Based on our land cover classification dendritic patterns of thermokarst gullies have formed throughout the valleys of the watershed. We estimate approximately 1/3 of the original 3.7 km 2 tundra la nd cover class has undergone near surface permafrost thaw that resulted in thermokarst formation. It is important to note that this is likely an underestimate of total permafrost degradation since thawing that does not result in subsidence (such as the deep en ing of active layer s) may not be remotely detected Thermokarst and tussock tundra classes occur in areas with relatively thick soil organic layers compared to other land cover cla sses within the EML watershed (F igure 2 3). Although our soil data is semi quantitative where we did not quantify variability >20 cm previous intensive work within tundra at EML reports soil organic layers range from 37 54 cm ( Hicks Pries et al. 2012 ; si tes occur within the subset of F igure 2 2). At these same sites, carbon inventories of the top meter of soil ranged from 55 to 69 kg C m 2 which is on the upper end of the range (16 94 kg C m 2 ) reported for similar tundra soils across Alaska ( Michaelson et al. 1996 ) They also found no difference in carbon pool s at minimally and extensively subsided tundra sites ( Hicks Pries et al. 2012 ) By extrapolating this intensive plot data, we estimate that both tundra and thermokarst classes have organic soil layers that are twice as deep as other classes in the watershed. In spite of the thermal protection provided by the extensive organic soil a ccumulations in the valleys ( Shur & Jorgenson 2007 ) under present climate conditions, this
26 protection is no longer sufficient to keep the underlying near sur face permafrost from thawing and thermokarst from forming and expanding The soil profile of thermokarst classes is seasonally thawing to a greater extent than any other land cover class, and thaw depth in thermokarst is twice as deep as in t he surrounding tussock tundra (F igure 2 4). This increase in thaw depth results from convection as water is redistributed into subsided areas as well as increases in soil thermal conductivity as organic matter becomes saturated. Both of these processes increase ground h eat flux within depressions while decreasing heat flux in higher, dryer areas ( Jorgenson et al. 2001 ; Kane et al. 2001 ; Bonan 2002 ; Yoshikawa & Hinzman 2003 ; Jorgenson & Osterkamp 2005 ; Osterkamp 2005 ; Osterkamp et al. 2009 ) The net result is deeper thawing in subsided, saturated areas and a decrease in thaw in the high relief areas of the adjacent tundra from which water is drained This pattern of increased thaw depth within subsided areas has been widely documented at EML ( Vogel et al. 2009 ; Lee et al. 2010 ; Belshe et al. 2012 ; Trucco et al. 2012 ) The landscape level changes created by permafrost thaw within the EML watershed have important implications for its carbon cycle because thermokarst features are forming in the valleys soils rich in organic carbon and are alter ing the physical environment in ways that increase the vulnerability of the soil carbon pool. Through associated changes in temperature and hydrology, thermokarst cause s the active layer of soil to thaw to a greater extent. This exposes more previously frozen organic carbon to above freezing temperatures and extends the time it remains unfrozen because it takes longer to refreeze the greater volume of unfrozen soil. These abiotic changes can stimulate decomposition and nitrogen mineralization, and result in increased carbon emissions, unless buffered by changes in soil hydrology ( Shaver et al. 1992 ; Goulden et al. 1998 ; Davidson & Janssens 2006 ; Schuur et al. 2008 ) On the other hand, thermokarst
27 formation can lead to increased carbon uptake if decompo sition increases plant available nitrogen ( Chapin & Shaver G. 1996 ; Hobbie et al. 2000 ; Ma ck et al. 2004 ) or the plant community shifts to more highly productive species ( Schu ur et al. 2007 ; Osterkamp et al. 2009 ) Currently, tundra undergoing permaf rost thaw at EML is a carbon source on an annual basis, emitting on average 54 g C m 2 year 1 over the past few years (2008 2010; ( Belshe et al. 2012 ) Al though both growing season net carbon uptake and winter carbon emissions have increase d in conjunction with ground subsidence at EML ( Vogel et al. 2009 ; Trucco et al. 2012 ) growing season net sink activity is negated by winter carbon emissions ( Be lshe et al. 2012 ) We believe this current imbalance is created because decomposition can continue to increase as permafrost thaw exposes more of the vast soil organic pool, while the magnitude of potential carbon uptake is mediated by size limitations of species currently residing in the tundra. This trajectory will continue as more soil is exposed with thaw, unless offset in the future, on the timescale of plant succession, by boreal trees advancing from the mixed riparian land cover class into the tun dra. An important next step is to build upon this analysis to detect changes in thermokarst extent as well as shifts in vegetation cover with remote sensing image time series Concluding Remarks Using a supervised maximum likelihood classification approach and a very high resolution IKONOS image, we detected and mapped small, irregular thermokarst features occur ring within the EML watershed (F igure 2 2). Overall we achieved an acceptable classification accuracy of 83%, with thermokarst classes having compar able accuracy to o ther classes in the landscape (T able 2 3). Our analysis shows that twelve percent of the landscape at EML has undergone subsidence and thermokarst induced physical changes have important implications for the carbon balance of this ecosys tem. The considerable time and work necessary to collect ground truth data within our relatively small (10.4 km 2 ) study area indicate that
28 mapping of such thermokarst features remains challenging despite their importance for understanding and projecting tundra ecosystem transitions. Homogeneous very high resolution remote sensing data will be useful to investigate such features in other areas. However, we believe that only a combination of very high resolution imagery and elevation data, such as from LIDAR, will allow successful mapping of such small, but widespread upland thermokarst features over larger study regions With the horizontal resolution (15 30m) and vertic al accuracy (>5 10 m) of most digital elevation models currently available in high latitude areas, it is not possible to extend this type of analysis to larger areas. Consequently, advances in other spatial products are needed, in conjunction with high res olution images, to detect formation of small thermokarst features that form due to climate or disturbance induced changes in high latitude ecosystems.
29 Table 2 1. Characterization of land cover classes and the number of field training sites u sed for the supervised classification. Land Cover Class Training sites Description Hilltop Shrub (HS) 10 Hill top communities dominated (85%) by medium stature (60 120 cm) shrubs ( Betula nana, Betula glandulosa ). *Similar to DLMP land cover class closed low shrub birch. §Similar to the CAVM land cover class S2 low shrub tundra. Shrub Tundra Transition (ST) 10 Intermediate transition zone between hilltops and valleys with intermediate stature (25 55 cm) vegetation that is a mixture of shrubs (63%; Betula nana Rubus chamaemorus Rhododendron groenlandicum Empetrum nigrum Vaccinium vitis idaea Vaccinium ulgino sum ), graminoids (16%; Eriophorum vaginatum Carex bigelowii ) and non vascular (17%; Dicranum spp., Sphagnum spp. Feather moss, Lichen) plants. §Similar to the CAVM land cover class S2 low shrub tundra and G4 tussock sedge, dwarf shrub, moss tundra. Tussock Tundra (TT) 10 Mixture of short stature (>25 cm) graminoid (35%; Eriophorum vaginatum Carex bigelowii ), shrub (38%; Betula nana Rubus chamaemorus Rhododendron groenlandicum Empetrum nigrum Vaccinium vitis idaea Vaccinium ulginosum ) and non va scular (25%; Dicranum spp., Sphagnum spp. Feather moss, Lichen) vegetation occurring in gently sloping valleys. *Similar to DLMP land cover class low shrub sedge. §Similar to the CAVM land cover class G4 tussock sedge, dwarf shrub, moss tundra. Thermoka rst Pool (TP) 4 Pools of standing water created by permafrost degradation in areas of low slope. Pools are generally around 5 m wide and less than 50 cm deep. Thermokarst (TK) 22 Thermokarst gullies that have three types of vegetation depending on the local slope and drainage. Relatively low slope areas are dominated by either non vascular species (60%; Sphagnum spp. ) or graminoid species (67%; Eriophorum spp ,Carex bigelowii ), while area wit h greater slope and drainage are dominated by a mixture of intermediate stature (20 63 cm) shrub species (61%; Betula nana Rubus chamaemorus, Rhododendron groenlandicum Empetrum nigrum Vaccinium vitis idaea Vaccinium ulginosum ). §Similar to the CAVM la nd cover class S2 low shrub tundra in high slope areas and W4 sedge, moss, low shrub wetland in low slope areas. Mixed Riparian (MR) 18 Mixed vegetation near the inlet and outlet of the lake, and alongside the road. Vegetation is a matrix of patches of needle leaf trees ( Picea mariana ) and large shrubs (200 600cm; either Salix spp or Alnus spp .). *Similar to DLMP land cover classes Stunted Spruce, Alder, and Willow. §Similar to the CAVM land cover class S2 low shrub tundra. Road (RD) 4 Stampede Road *Similar to DLMP land cover class bare ground. Lake (LA) NA Eight Mile Lake. *Similar to DLMP land cover class clear water. Denali Landcover Mapping Project (DLMP); Stephens et al. 2001 § Circumpolar Arctic vegetation map (CAVM); Walker et al. 2005
30 Table 2 2. Spatial extent and percent of total area covered by each class within the study area. Land Cover Class Area % Area km 2 Hilltop Shrub (HS) 15.8 1.7 Shrub Tundra Transition (ST) 22.4 2.3 Tussock Tundra (TT) 22.9 2.4 Thermokarst (TK) 12.0 1.3 Thermokarst Pool (TP) 0.1 0.01 Mixed Riparian (MR) 21.0 2.2 Lake (LA) 5.2 0.5 Road (RD) 0.5 0.1 Total Area 100 10.4
31 Table 2 3 Summary of supervised classification accuracy assessment that includes a confusion matrix of the classification and ground truth data, along with the producer's and user's accuracies (%), the overall accuracy, and the Kappa statistic. Ground truth (pixel) Land Cover HS ST TT TK TP MR LA RD Total Producer's User's HS 1501 234 0 61 0 166 0 0 1962 85 74 ST 60 3999 216 100 0 7 0 0 4382 69 77 TT 1 126 8624 7 0 0 0 0 8758 97 88 TK 78 538 66 459 0 23 0 0 1164 79 75 TP 0 0 0 2 45 2 2 0 51 83 100 MR 68 4 0 19 0 338 0 0 429 77 90 LA 0 0 0 0 0 0 836 0 836 100 100 RD 0 0 0 0 0 0 5 82 87 100 100 Total 1708 4901 8906 648 45 536 843 82 17669 Overall Accuracy 83% Overall Kappa 0.79
32 Figure 2 1. IKONOS image taken during June 2008 of the EML watershed and adjacent hill slope (10.4 km 2 study area outlined in white) located in Healy, Alaska. Training sites used for the supervised classification and validation sites used for accuracy assessment are a lso shown.
33 Figure 2 2. Mapped supervised classification of the EML watershed and adjacent hill slope, with a boxed subset showing the area where previous studies measured carbon pools and fluxes.
34 Figure 2 3. Soil organic layer thickness (cm) of l and cover classes, with colors corresponding to the mapped supervised classification. Classes with organic depths less than 20 cm (MR, HS, ST) are shown in box plots, with solid line denoting mean and box denoted the 25 th & 75 th quartile. Classes grouped i nto the >20cm category (TT, TK, TP) have boxes that start at 20cm and fade toward 40 cm to represent that previous site measurements show organic depths are greater than 20cm (ranging from 37 54 cm) for these classes.
35 Figure 2 4. Maximum seasonal thaw depth (cm) of land cover classes, with solid line denoting mean and box denoted the 25 th & 75 th quartile, and colors corresponding to the mapped supervised classification. Letters denote differences among classes at a
36 CHAPTER 3 INCORPORATING SPATIAL HETEROGENEITY CREATED BY PERMAFROST THAW INTO A LANDSCAPE CARBON ESTIMATE Background Northern high latitudes are disproportionally warming and Arctic temperatures are predicted to increase by 6.5 C or more by the year 2100 in response to radiative forcing caused by increasing greenhouse gases and changes in albedo ( Chapin et al. 2000 ; Chapin et al. 2005 ; Hinzman et al. 2005 ; IPCC 2007 ) Currently, permafrost occur s within 24% of the ice free land area in the northern hemisphere ( Zhang et al. 1999 ) and it is estimated 25 90% will degrade into seasonally frozen ground by the year 2100 ( Anisimov & Nelson 1996 ; Saito et al. 2007 ; Lawrence et al. 2008 ) According to recent estimates, permafrost soils contain twice as much carb on (1672 Pg) as the entire atmospheric pool ( Schuur et al. 2008 ; Tarnocai et al. 2009 ) If a portion of this C is released to the atmosphere it could result in a strong positive feedback to climate change. Understanding how permafrost thaw affects the rate of C exchange from this large pool is essential for understanding the global C cycle in a warmer world. Thawing of permafrost is a temporally dynamic and spatially heterogeneous process. Rising temperatures increase active layer (seasonally thawed surface layer) thickness and form thermokarst ( Jorgenson & Osterkamp 2005 ; Zhang et al. 2005 ) Thermokarst is uneven ground that forms when ice rich permafrost thaws, drainage occurs, and the ground surface subsides ( Jorgenson & Osterkamp 2005 ) These localized changes in surface relief greatly alter the surface hydrology of the area. As water is redistributed from higher to lowe r microtopographical areas, thermal erosion by the movement of water warms the soil and further perpetuates Belshe EF, BM Bolker, R Bracho, EAG Schuur (2012) Incorporating Spatial Heterogeneity Created by Permafrost Thaw into a Landscape Carbon estimate. Journal of Geophysical Research: Biogeosciences Vol. 117.
37 permafrost thawing ( Kane et al. 2001 ; Osterkamp et al. 2009 ) This positive feedback creates a mosaic of patches that range from high, dry embankments with shallow active layers to subsided areas with relatively wet, warm soils and deep active layers ( Osterkamp et al. 2009 ; Vogel et al. 2009 ; Lee et al. 2010 ) Furthermore, this pattern of ground subsidence, dictated by the initial p resence of ice rich permafrost, is interspersed throughout the landscape and ultimately creates a mosaic of various degrees of permafrost thaw and microtopography The future C balance of high latitude ecosystems depends on the sensitivity of biological pr ocesses (photosynthesis and respiration) to the physical changes in temperature and moisture occurring with permafrost thaw. But, predicting C exchange in these ecosystems is difficult because of the landscape heterogeneity created as permafrost thaws. Sin ce adjacent patches can have very different physical environments, they can have very different gross primary production (GPP) and ecosystem respiration (R eco ) ( Vogel et al. 2009 ; Lee et al. 2010 ) Landscape scale GPP and R eco will depend on the cumulative response of the landscape to permafrost thaw, which in turn will dictate the direction and magnitude of net ecosystem exchange (NEE = GPP R eco ). The respon se of the C cycle to spatial and temporal environmental variation is often non linear and not simply described by the mean response ( Aubinet et al. 2002 ) Therefore, an appropriate understanding of both the spatial and temporal variation of C flux is essential for estimating the C balance of a landscape. But the intensive temporal sampling required for good estimates of C flux makes it difficult to obtain extensive, sp atially explicit C flux data. Eddy covariance (EC) provides a method to directly measure C exchange at a high level of temporal resolution over a large spatial scale ( Baldocchi 2003 ) But fluxes measured by EC are commonly assumed to come from a homogenous surface, which makes it difficult to resolve the cumulative contribution of localized features in the landscape to an EC estimate ( Schmid & Lloyd 1998 ;
38 Laine et al. 2006 ) Although much effort has gone into developing and using models to locate where fluxes originate (i.e. the footprint of an EC tower) ( Schmid 1997 ; Schmid & Lloyd 1998 ; Kormann & Meixner 2001 ; Schmid 2002 ) less effort has gone into incorporating spatial information b ack into EC C estimates. However, the abundance of data produced by EC towers gives us the ability to explore spatial patterns of C flux. In this study, we use generalized additive models (GAMs) to generate a continuous time series of NEE for a tundra land scape undergoing permafrost thaw. We developed a spatially explicit quantitative metric of permafrost thaw based on variation in microtopography. By incorporating our spatial metric into EC gap filling models, we were able to make C flux predictions for th e landscape as a whole, as well as for specific landscape patches throughout the continuum of permafrost thaw and ground subsidence. We tested the robustness of our models against more widely used (non spatial) gap filling methods. Our objectives were to m ore accurately estimate the C balance of a heterogeneous landscape and to explore the collective effect of permafrost thaw on the plant and soil processes that dictate ecosystem C exchange. Material and Methods Site Description the northern foothills of the Alaska Range near Denali National Park and Preserve ( Schuur et al. 2007 ; Schuur et al. 2009 ) This upland area occurs within a vulnerable band of discontinuous permafrost near the point of tha w due to the combination of its elevation and geographic position ( Yocum et al. 2006 ; Romanovsky et al. 2007 ) Deep permafrost temperature has been measured at the site since 1985 and during this time thermokarst terrain has developed and expanded as the permafrost has warmed ( Osterkamp et al. 2009 ) Vegetation at the site is dominated by moist acidic tussock tundra comprised of sedge ( Eriophorum vaginatum), deciduous and evergreen
39 shrubs ( Vaccinium uliginosum, Rubus chamaemorus, Betula nana and Ledum palustre ), and non vascular plants ( Sp hagnum spp ., Dicranum spp., feathermoss and lichens). Soils at the site are classified as Gelisols because permafrost is found within 1 meter of the soil surface ( Soil Survey Staff 1999 ) An organic horizon, 0.45 0.65 m thick, covers cryoturbated mineral soil that is a mixture of glacial till (small stones and cobbles) and windblown loess. Organic C pools in the top meter of soil range between 55 69 kg C m 2 ( Hicks Pries et al. 2012 ) The long term mean annual air temperature (1976 2009) of the area is 1.0 C and the g rowing season (May Sept) mean air temperature is 11.2 C, with monthly averages ranging from 16 C in December to +15 C in July. The long term annual mean precipitation is 378 mm with a growing season mean precipitation of 245 mm (National Climatic Data Center, NOAA) Mean growing season air temperature was 8.1 C and 9.7 C during 2008 and 2009, respectively, and growing season precipitation was 346 mm and 178 mm during 2008 and 2009, respectively. Eddy Covariance Measurements NEE was measured using eddy covariance (EC) from June to December 2008 and April to October 2009. The EC system consisted of a CSAT3 sonic anemometer (Campbell Scientific, Logan, UT) and an open path CO 2 /H 2 0 gas analyzer (Li 7500, LI COR Biosciences, Lincoln, NE) mounted on a 2m tow er. Data were recorded at a frequency of 10 Hz on a CR5000 data logger (Campbell Scientific, Logan, UT), and fluxes were Reynolds averaged over 30 min time periods ( Reynolds 1985 ) Calibration was preformed monthly during the growing season using a zero CO 2 air source, a 1% standard CO 2 concentration, and a dew point generator (Li 610, LI COR Biosciences, Lincoln, NE) for water vapor. The EC tower was pla ced within a patchy landscape consisting of visibly subsided areas to the North and West and relatively even terrain to the South and East. The fetch from the tower was greater than 300 m in all directions and winds predominantly came from the NE and SW. A n analytical footprint model developed by
40 Korman and Meinxner (2001) showed on average 50% of fluxes originated within the first 50 m around the tower, and greater than 80% of fluxes originated within 200 m from the tower. Eddy Covariance Data Handling Raw CO 2 fluxes were corrected for damping of high frequency fluctuations, sensor separation, and misalignment of wind sensors with respect to the local streamline ( Moncrieff et al. 1997 ; Aubinet et al. 2000 ; Wilczak et al. 2001 ) CO 2 fluxes were then corrected for variations in air density due to fluctuation in water vapor and heat fluxes ( Webb et al. 1980 ) and for fluctuations caused by surface heat exchange from the open path sensor during wintertime conditions ( Burba et al. 2008 ) Data screening was applied to eliminate half hourly fluxes with systematic errors and non relevant environmental influences such as: (1) incomplete half hour data sets as a result of system calibration or maintenance, (2) time periods when the canopy was poorly coupled with the external atmospheric conditions as defined by the friction velocity, u (threshold < 0.12 m s 1 ; ( Goulden et al. 1996 ; Clark et al. 1999 ) (3) excessive variation from the half hourly mean based on an analysis of standard deviations for u, v, and w wind statistics and CO 2 fluxes. Fluxes were then divided into weekly data sets for both day a nd night conditions and unrealistic low or high values (>2 standard deviations fro m the mean) were filtered out. In total, ecosystem fluxes were measured 72% and 96% of the time during the 2008 and 2009 campaigns, respectively, while 64% and 60% of those v alues were eliminated by the screening criteria listed above. The quality of our data was evaluated by the degree of growing season energy closure (Rnet = LE + H + G), which was 76% in 2008 and 73% in 2009. Ground heat flux (G) was estimated as the change in soil temperature with depth plus soil heat storage ( Liebethal et al. 2005 ; Liebethal & Foken 2007 ) To calculate soil heat storage, we assumed 40% organic matter content ( Hicks Pries et al. 2012 ) and 60% volumetric water content based on soil cores taken from the site. Measurements of half hour NEE were calculated as: NEE = Fco 2 + Fs, where Fco 2
41 was the mean flux of CO 2 at measurement height and Fs was the half hour change in CO 2 stored below meas urement height. Because of the short vegetation (~30 cm), we calculated the change in CO 2 storage by taking the difference in successive CO 2 measurements at the measurement height ( Hollinger et al. 1994 ) We used the meteorological convention that positive NEE represents a transfer of CO 2 from the ecosystem to the atmosphere. Environ mental Measurements Standard meteorological data were collected on a tower adjacent to the EC tower, including photosynthetic photon flux density (PPFD; Li 190SA, Li Cor Inc. Lincoln, NB.), incident radiation (Li 200SA, Li Cor Inc. Lincoln, NB), net radiat ion (REBS Q*7.1, REBS. Inc., Seattle, WA.), relative humidity and air temperature (Vaisala HMP45c, Campbell Scientific Inc., Logan, UT), and wind speed and direction (RM Young 3001, Campbel l Scientific Inc., Logan, UT). Soil temperature profiles (5, 10, 15 20 and 25 cm from surface) were measured with constantan copper thermocouples and a thermistor (at 5 cm depth only; 107, Ca mpbell Scientific, Logan, UT). Moisture integrated over the top ~15cm of soil was measured with a Campbell CS615 water content refl ectometer. All measurements were recorded at half hour average intervals with a CR5000 data logger (Campbell Scientific, Logan, UT). A complete replicate set of micrometerological measurements were collected at a tower 100 m to the NW of EC tower, and were used to interpolate gaps in micrometeorological data measured at the EC tower. Landscape Properties To quantify the amount and distribution of land surface subsidence associated with permafrost thaw, a digital elevation model (DEM) was created from point measurements of elevation. Fine scale differences in elevation were measured with a high resolution differential global position system (dGPS). One GPS unit (Trimble 5400) was placed at a nearby USGS which acted as the reference
42 receiver. Using a second GPS unit (Trimble 5400) secured to a backpack, a kinematic survey was conducted by walking transects within a 400 m diameter circle encompassing the EC tower footprint. Geographic position and elevation were collected at 5 s intervals, yielding a total of 7,220 points. These data were post processed with methodology developed by UNAVCO using Trimble Geomatics Office (Dayton, Ohio). To create the DEM of the area surrounding the EC tower, spherical models were fit to empirical semivariograms, and ordinary kriging was used to interpolate between point measurements using the calculated range of 282 m, a nugget of 0.02, and a partial sill of 0.47. Variogram analysis and kriging was done with t he Geostatistical analyst extension in ArcGIS 9.3. Because the study site is on a gentle slope (~5%), the original DEM was corrected for overall slope elevation changes, so we could decipher small scale subsidence features. To correct for the slope, the DE M was first rescaled so the minimum elevation equaled zero. Mean elevation within 30 m blocks was subtracted from the rescaled DEM, resulting in the deviation in elevation away from the mean plane. This created a map of small scale variations in topography that we d efine here as microtopography. Pixel resampling and calculations were done using the aggregate function, resample tool, and the raster calculator in the Spatial Analyst extension in ArcGIS 9.3. To obtain landscape information in a form comparable to EC data, we extracted information on microtopography corresponding to each wind direction sampled by the EC tower. Virtual transects 200 m in length, originating at the EC tower and radiating out in every wind direction (0 359), were created. A distanc e of 200 m was chosen because it corresponded to the distance where on average >80% of scalar fluxes originated based on an analytical footprint model ( Kormann & Meixner 2001 ) Microtopography (i.e., local elevation) was sampled every
43 ( Beyer 2004 ) in ArcGIS 9.3. The standard deviation of microtopography, which we refer to as roughness, was calculated for each transect (wind direction). This calculated metric was chosen because it captures the variation in microtopography created by permafrost thaw (both subsided areas and raised embankments). Our metric, roughness, should not be confused with the micrometeorological term roughness length. To calculate our metric, roughness, corresponding to each (half hour) flux measurement, we simply calculated the standard deviation of the per meter values of microtopography along the entire transect corresponding to the measured wind direction. We acknowledge that C fluxes measured over a 30 minute period do not emanate f rom a one dimensional transect; instead they come from two dimensional areas in the landscape. To find the best spatial metric corresponding to the measured C flux, we also calculated roughness for the 3 and 5 adjacent transects of the measured wind direct ion. We found no change in the relationship between C flux and roughness when using a greater number of transects. Also, in principle, footprint models provide more information than the overall radial scale of the area surrounding the EC tower because they help to pinpoint where in the landscape fluxes are originating ( Schmid 1997 ; Schmid & Lloyd 1998 ; Kormann & Meixner 2001 ; Schmid 2002 ) So, for comparison, we also used estimates of the cumulative probability of fluxes coming from different fetches, calculated by a footprint model, to calculate a weighted standard deviation of roughness. However, we chose to use the simple non weighted roughness because under certain conditions weighting caused relatively flat areas to the SE to have a higher standard deviation than the most subsided areas to the NW. We believe thi s discrepancy is due to the mismatch in scale between our one dimensional transects and the two dimensional cumulative density function calculated by the footprint model ( Kormann & Meixner 2001 ; Kljun et al. 2003 )
44 We explored the relationship between roughness and normalized difference vegetation index (NDVI), and the relationship between microtopograghy and active layer depth (ALD). N DVI was calculated using spectral data from an IKONOS image of the site acquired in June 2008. Mean NDVI was calculated for each of the 360 virtual transects radiating out from the EC tower, and was compared to the roughness of the corresponding transect. ALD was measured at 310 locations stratified at various distances within the potential EC footprint, by measuring the length of a metal probe inserted into the soil until the impenetrable frozen layer was reached. The geographic location of each site was m easured and subsequently used to extract corresponding values of elevation from the map of microtopography. We only compared ALD to microtopography because we did not have a continuous surface of ALD; therefore, were unable to extract data for the 360 virt ual transects to compar e with roughness. Relationships were explored with generalized additive models using the mgcv package in R ( Wood 2008 ; R Development Core Team 2010 ) Estimation of Landscape Scale Carbon Exchange To estimate the carbon balance of an ecosystem, missing values in measured CO 2 flux time series must be gap filled to generate a continuous time series of net ecosystem exchange (NEE). We estimated carbon e xchange using two gap filling strategies: 1) a novel gap filling strategy using generalized additive models (GAMs) that are flexible enough to incorporate spatial information; 2) nonlinear (NL) relationships with non spatial environmental variables. Althou gh NEE is directly measured by the EC technique, the driving force of the exchange is dependent on environmental conditions. Therefore, we modeled NEE for gap filling during winter, growing season (GS) days, and GS nights separately. The beginning and end of the GS was determined by abrupt changes in net radiation corresponding to snowmelt and widespread snow cover, respectively. Generally, the GS began in early May and ended at the end
45 of September. Data during the GS were split into day and night by ambie nt light, so when PPFD was greater than 10 mol m 2 s 1 daytime conditions were assumed. During daytime, NEE is the balance between gross primary production (GPP) and ecosystem respiration (R eco ). To tease apart their contributions, we modeled R eco during GS days using models fitted with GS night data and calculated GPP as the difference between NEE and R eco (GPP = NEE R eco ). Once soil temperature at 5 cm fell below 0 C, winter conditions were assumed. During these conditions, photosynthesis is no t occurring so NEE is equivalent to R eco Gap filling strategy 1: generalized additive models To generate a continuous time series, we gap filled NEE using generalized additive models (GAMs), an extension of generalized linear models where a response is mo deled as the additive sum of smoothed covariate functions ( Hastie & Tibshirani 1990 ; Wood 2006 ) With GAMs, nonlinear effects can be modeled without manually specifying the shape of the relationships, which provided us the flexibility to incorporate roughness along w ith other explanatory variables into the prediction of NEE ( Wood 2006 ; Zuur et al. 2009 ) To contr ol the shape of functions, we used penalized regression splines, which determine the appropriate degree of smoothness of each smoothing function by generalized cross validation (GCV) and adds a smooth with maximum likelihood ( Wood 2006 ) All GAMs used had the basic form: where y i denotes the response variable (NEE or R eco 0 the intercept, and functions f 1 (x i ) and f 2 (z i ) are smooth functions of explanatory variables x i and z i and f 3 (x i z i ) is a two dimensional smooth function of their interactions. We used thin plate regression splines as the basis for representing smooths ( f 1 and f 2 ) for single covariates and tensor pr oduct smooths ( f 3 ) for
46 interactions (multiple covariates) because they have been found to perform better when covariates are not on the same scale ( Wood 2006 ) We forced the effective degrees of freedom in each model to count as 1.4 degrees of freedom in the GCV score, which forces the model to be slightly smoother than it might otherwise be, this is an ad hoc way to avoiding overfitting ( Kim & Gu 2004 ) Because eddy covariance data are heteroscedastic ( Richardson et al. 2008 ) and there were distinct patterns in the residuals, we fit GAMs using a mixed model framework (GAMM) with a Gaussian error distribution to facilitate the incorporation of an exponential variance structure: where the variance of the r 2 is multiplied by an exponential function of the fitted values ( V) We used all data in a single dummy group as our random effect to facilitate the incorporation of the variance structure into the GAM ( Wood 2006 ; Dormann 2007 ; Zuur et al. 2009 ) All models were fitted using the mgcv package ( Wood 2006 ) in R ( R Development Core Team 2010 ) A subset of explanat ory variables was selected a priori including: PPFD, temperature (air, soil at 5cm, depth integrated soil temperature down to 25cm), roughness, and day of year (DOY). We suspected there might be complex interactions between explanatory variables, so both d irect effects and all possible interactions were compared. Models were selected for each time period (Winter, GS day, GS night) during each year (2008, 2009), by starting with the full model containing all variables and interactions and using a form of aut omatic backward selection in which the penalization term for each smooth could automatically set the term to zero and remove it from the model as appropriate ( Wood 2008 ) We also took into consideration how removing terms affected: (1) the GCV score (the lower the better); (2) the deviance explained
47 (the higher the better) and (3) the Akaike Information Criterion (AIC, the lower the better; ( Anderson et al. 1998 2001 ) Because our GAM models incorporated landscape information, they allowed us to estimate the landscapes carbon ba lance in two different ways. If we assumed the land scape was filled depending on the measured wind direction at the time of the gap. This resulted in a single time series of carbon exchange for the landscape (G AM 1). Alternately, if we assumed the landscape was a combination of multiple patches (wind directions), all absorbing or releasing C simultaneously, then each wind direction was gap filled separately for the entire time series. This resulted in 360 separa te time series, whose predictions were averaged to achieve an estimate of carbon exchange for the entire landscape (GAM 360). This method allowed us to estimate C exchange for the entire heterogeneous landscape and made it possible to compare predictions f rom landscape patches that differed in roughness. Both methods of prediction were done for each time period during 2008 and 2009. Gap filling strategy 2: non linear regressions For comparison we also gap filled data using a more traditional non li near (NL) regression approach. During GS days, gaps were filled using parameters obtained by fitting half hour NEE to PPFD using a non rectangular hyperbola ( Thornley & Johnson 1990 ) : active radiation, P max is the asymptote, and R is the intercept or dark respiration term. To capture changes in phenology, parameters were estimated biweekly or monthly depending on the variation among weeks. We incorporated an exponential variance struct ure due to the heteroscedacity of the data
48 and used maximum likelihood to estimate parameters using the bbmle package ( Bolker 2010 ) in R. We were unable to fit exponential models to winter and GS night data separately, so we gap filled with parameters estimated using both data together. Parameters were estimated using the following equation: integrated soil temperature during 2008 and soil temperature measured at 5 cm during 2009. We compared models with various forms of temperature (air, soil at 5 cm, depth integrated soil) and chose the best model based on than the alternatives. An exponential variance structure was added and parameters were estimated using generalized non linear least squares using the nlme package in R ( Pinheiro et al. 2011 ) Model Performance We compared the coefficient of variation (R 2 ) and Akaike Information Criterion (AIC, ( Anderson et al. 1998 2001 ) of each GAM and NL model, during each time period. Because of the variation in model types, w e calculated the R 2 simply as the correlation between the predicted values from each model and the observed values. To assess the predictive performance of the GAM and NL models, we performed cross validation. Ten percent of the data was randomly removed, models were fitted to the remaining data, and these models were then used to predict responses for the withdrawn ten percent. This process was repeated ten times and the root mean square error (RMSE) was calculated for each model. We then compared RMSE of the GAM and NL models within each time period using a t test at a statistical significance of P<0.05. To calculate comparable values of AIC we used the following equation:
49 where n=number of observations and p=number of parameters ( Venables & Ripley 2002 ) Results Landscape Heterogeneity of the EC Footprint Our map of microtopography captured the spatial pattern of ground subsidence created by permafrost tha w within the EML watershed (Figure 2 1). The largest variation in microtopography was found to the NW of the E C tower, while areas to the E and SE were relatively flat. This spatial distribution of ground subsidence agreed well with patterns visible in high resolution aerial photographs of the site. In addition, the maximum (0.5 m) and minimum ( 0.9 m) deviations away from the mean elevation were consistent with field measurements of the depth of individual subsided features and height of raised embankments created by thaw (data not shown). This pattern was mirrored by our calculated landscape metric roughness. Tra nsects with the highest and lowest roughness corresponded with the winds coming from the N NW and SW SE, respectively (Figure 3 1). In general, the depth of the active layer increased as local elevation decreased (became more subsided) and microtopography explained 51% (adj. R 2 =0.51) of the observed variation in active layer depth. The relationship was non linear, with little variation in active layer depth at sites where elevation was positive or slightly negative, followed by an exponential increase in ac tive layer depth as e leva tion fell below 0.2 m (Figure 3 2A). Transects with more variation in microtopography (high roughness) were found to have higher mean NDVI than transects with less variation (low roughness), and roughness explained 55% of the vari ation in mean NDVI (adj. R 2 =0.55, edf=6.1). This relationship was also non linear, with NDVI linearly increasing with roughness, then leveling out and sometimes decreasing as roughness increased a bove 0.14 (Figure 3 2B).
50 The magnitude of net ecosystem exch ange (NEE) increased as the roughness of the landscape increased (Figure 3 3). During GS days, NEE became more negative (more C uptake) with increasing roughness. More C was released from areas with higher roughness during GS nights, but the trend of incre ased C emission with increasing roughness decreased in magnitude and then reverse d during the winter months (Figure 3 3). Model Performance GAM models outperformed the non spatial NL models for gap filling C exchange. GAM models had a higher or equivalent coefficient of determination (R 2 ) and higher predictive power (lower AIC) than NL models during every time period in both 2008 and 2009 (Table 3 1). During cross validation, GAM models always had lower mean RMSE, but were only significantly lower (at p<0.0 5) during GS days of 2008 and 2009 (Table 3 1). Predictions of Ecosystem C Balance In general, the three gap filling methods (NL, GAM 1, GAM 360) resulted in similar estimates of NEE, GPP and R eco for the time periods June 6 through December 8 2008 (weeks 24 48) and April 24 through Octob er 10 of 2009 (weeks 12 40; Figure 3 4). Weekly estimates of NEE, GPP and R eco generated for both GAM methods closely mirror one another thr oughout both 2008 and 2009 (Figure 3 4). Although the NL and two GAM methods generated similar final estimates of net ecosystem exchange, there were some notable differences. The two GAM methods estimated a slightly higher uptake of carbon (more negative GPP) than their NL counterpart, 10 13 gCm 2 more in 2008 and 12 19 gCm 2 m ore in 2009. This difference in GPP was spread throughout the growing season, with no single week solely resp onsible for the difference (Figure 3 4A). Similarly, the two GAM methods estimated a slightly higher release of carbon (R eco ) than the NL method, 1 0 11 gCm 2 more in 2008 and 7 13 gCm 2 more during 2009. Unlike GPP, however, differences in R eco could be attributed to certain time periods.
51 During 2008, the major differences between the predictions of the two methods occurred during the early growin g season (weeks 21 29) where the NL method predicted lower R eco than either GAM method. During the majority of the GS of 2009 the GAM methods estimated higher R eco than the NL method. The other time of divergence among the methods was during transitions in to and out of the growing season. The GAM methods predicted lower R eco than the NL method during transition from GS to winter (weeks 39 40) in 2008 and during the transition from winter to GS ( week 15 19) in 2009 (Figure 3 4B). Predictions of Landscape Het erogeneity of C Flux Using GAMs, we were also able to predict C exchange for each wind direction. To estimate the C balance for the entire landscape we calculated the mean carbon flux of all wind directions for each 30 min interval throughout 2008 and 2009 (GAM 360). This resulted in an estimate of NEE, GPP and R eco for the landscape on average but also allowed us to compare C fluxes from the wind directions with the minimum and maximum roughness to further understand the influence of permafrost thaw and gr ound subsidence on C flux. From June to December 2008, GAM 360 estimated the landscape on average took up 337.1 gCm 2 via photosynthesis and released 289.5 gCm 2 via respiration, resulting in an ecosystem carbon gain of 47.5 gCm 2 (Figure 3 5 and Table 3 2 ). The direction with the maximum roughness had higher GPP and R eco than the landscape on average, while the direction with the minimum roughness had lower GPP and R eco (Figure 3 5 and Table 3 2). This resulted in the direction with maximum roughness gaini ng 55% more C than the landscape on average, while the direction with minimum roughness gained 76.4% less C (Table 3 2). From April to October 2009, the landscape on average took up 498.7 gCm 2 via photosynthesis and released 410.3 gCm 2 via respiration, r esulting in a net gain of 87.8 gCm 2 (Figure 3 5 and Table 3 2). Again, the direction with the maximum roughness had higher GPP
52 and R eco than the landscape on average, while the direction with the minimum roughness had lower GPP and R eco (Figure 3 5 and Table 3 2). This resulted in the direction with maximum roughness gaining 41.4% more C than the landscape on average, while the direction with minimum roughness gained 61.7% less C (Table 3 2). The amplified increase in GPP with maximum roughness was consistent throughout all weeks of the growing se ason in both 2008 and 2009 (Figure 3 5A). Unlike GPP, the increase in R eco with roughness was not consistent throughout either year. During the GS of both years, the landscape with maximum roughness had hig her R eco but during the winter this trend reversed and the landscape with min roughness had higher R eco (Figure 3 5B). Even though areas with minimum roughness had higher R eco during winter, the overall carbon emission throughout 2008 and 2009 was still g reater from areas with maximum roughness (Table 3 2). Discussion Quantifying Spatial Heterogeneity Microtopography is an easy to obtain, integrative metric of the physical and biological changes occurring as the result of permafrost thaw within the EML watershed because it correlates with v ariables that drive C cycling. Roughness, our landscape level metric of permafr ost thaw, captured the variation in microtopography of each wind directi on sampled by the EC tower (Figure 3 1). We found that as microtopography decreased (ground became more subsided), active layer depths (ALD) increased, exponentially in creasing after a threshold (Figure 3 2A). This pattern is a result of changes in soil thermal conductivity created by the redistribution of water into subsided areas, which increases soil temperature within depressions while decreasing temperatures in higher, dryer areas ( Jorgenson et al. 2001 ; Kane et al. 2001 ; Osterkamp et al. 2009 ) These physical changes in soil moisture and temperature drive va riable depths of thaw across the hillslope. This landscape level pattern of ALD and subsidence is
53 consistent with previous work at this site, which showed similar relationships between microtopography, temperature, and moisture ( Vogel et al. 2009 ; Lee et al. 2010 ) Areas with greater roughness had higher mean NDVI, with NDVI increasing with roughness unt il a threshold was reached (Figure 3 2B). Unlike ALD, NDVI leveled out and slightly decreased at the upper end of the roughness scale. NDVI has been shown to be positively correlated with leaf area index ( Tucker 1979 ; Williams et al. 2008 ) aboveground biomass ( Sellers 1987 ; Boelman et al. 2003 ) net primary production ( Goward et al. 1985 ) GPP and R eco ( Boelman et al. 2003 ; Vourlitis et al. 2003 ; La Puma et al. 2007 ) Permafrost thaw within the EML watershed causes a shift in species composition from a plant community dominated by tussock forming sedges to a community with increased shrub and moss abundance, and concurrently an increase in biomass and productivity ( Schuur et al. 2007 ; Vogel et al. 2009 ) Our result of increased NDVI with increased roughness is consistent with this pattern of increased biomass and productivity with thaw and also indicates that an upper limit of productivit y may be reached as permafrost thaw continues and plants respond to the changing conditions. This upper limit is likely driven by the size of shrub species currently at the site, but could increase in the future, on the timescale of plant succession, if bo real trees were to increase in abundance at this tundra site. These relationships between microtopography and important biophysical features of the landscape (ALD & NDVI) highlight the feasibility of using remotely sensed spatial information to improve est imates of regional C balance in high latitude ecosystems. Recent advancements in sensor resolution (e.g. LIDAR) now make microtopographical mapping of these vast, remote areas possible.
54 Incorporating Spatial Heterogeneity into EC C Estimates We incorporate d the spatial variability of C flux into the EC estimate of C exchange during gap filling using generalized additive models (GAMs). Because a continuous time series is required to estimate C balance, missing time periods must be modeled ( Falge et al. 2001 ; Baldocchi 2003 ) This gap filling step provided a method for predicting C exchange based on the roughness of the landscape in specific wind directions, as well as more accurat ely determining the C balance of the entire landscape. We found that GAMs were equivalent or superior to traditional NL regression approaches (Table 3 1). GAMs had higher predictive power and a higher or equivalent coefficient of variation (R 2 ) than the NL models during all time periods. During cross validation, GAMs consistently had lower RMSE than NL models over all time periods, but were statistically lower only during GS days. The lack of statistical improvement in RMSE during GS nights and winter by th e GAMs is not surprising because these time periods are notoriously difficult to model using any procedure ( Baldocc hi 2003 ) Overall, our model comparison to the more traditionally used NL gap filling models gave us confidence that adding model complexity to include spatial information was justified. The aggregated predictions of C exchange from the GAM and NL model s did not substantially differ from one another thro ughout either 2008 or 2009 (Figure 3 4). However, there were notable time periods where the two methods diverged in their predictive capabilities. During the early growing season of 2008 (weeks 21 29), the GAM substantially over predicted R eco compared to the NL model (Figure 3 4B). We believ e this difference is due to a lack of data during the early GS, which caused the GAM to miss the upswing of R eco that coincides with rapid changes in phenology during the spring. Predictions of R eco and NEE from the NL and GAM also diverged during the tran sition from the GS to winter in 2008, and the transition from the winte r to the GS of 2009 (Figure 3 4B & Figure 3 4C). We attributed this sensitivity to
55 seasonal transitions by the GAM to the flexibility of their smoothing functions, which can capture rap idly changing trends in the data ( Zuur et al. 2009 ) Because the empirically derived parameter estimates depend on relationships that change dramatically in this highly seasonal ecosystem, the relationships would need to be continuously updated in order to capture these transitions ( Falge et al. 2001 ; Baldocchi 2003 ) Biologically, the dip in R eco during the transition in and out of winter could be attributed to changes in microbial species composition, or to the disruption of biological activity by the state change (freezing point) of water ( Rivkina et al. 2000 ; Mikan et al. 2002 ) This dip in R eco could also be caused by shifts in the availability and use of substrates by microbes between the GS and winter ( Hobbie et al. 2000 ; Dioumaeva et al. 2003 ; Davidson & Janssens 2006 ) Effects of Spatial Heterogeneity on C Flux By incorporating spatial information, we were able to estimate the C balance of the landscape in two different ways. First, we filled gaps depending on the measured wind direction at the t ime of the gap and created a single time series of C exchange for the entire landscape (GAM 1). Second, we gap filled the entire time series for each wind direction separately and averaged the predictions to estimate the C balance of the landscape (GAM 360 ). This allowed us to estimate C exchange for the entire heterogeneous landscape and also make predictions for the wind directions with the minimum and maximum roughness. Final C estimates from GAM 1 and GAM 360 were nearly identical (Figure 3 4). This similarity indicates that the wind distribution sampled by the EC tower was sufficient to capture the variability of thaw seen across the entire radial landscape surrounding the tower. The landscape on average (GAM 360) was a C sink during both 6 month mea surement campaigns. The wind direction with the most variation in microtopography (maximum roughness), resulting from permafrost thaw and ground subsidence, had both higher GPP and
56 R eco tha n the landscape on average (Figure 3 5A & Figure 3 5B), while the w ind direction with least variation (minimum roughness) exhibited lower GPP and lower R eco Overall, during the 6 month campaign of 2008, the area with highest roughness gained 55.2% more C than the landscape on average, while the area with lowest roughness gained 76.4% less C. Similarly, during the 6 month campaign of 2009, the area with the highest roughness gained 41.4% more C, while the areas with lowest roughness gained 61.7% less C (Table 3 2). Based on these results, permafrost thaw and ground subside nce amplifies both GPP and R eco The amplification of GPP with roughness was consistent throughout th e GS of both 2008 and 2009 (Figure 3 5A & Figure 3 3). This enhanced C sequestration could be due to shifts in the plant community to more highly productiv e species as permafrost thaws ( Schuur et al. 20 07 ; Osterkamp et al. 2009 ) or to increased plant productivity due to greater nutrient availability resulting from enhanced decomposition within subsided areas ( Shaver et al. 1992 ; Mack et al. 2004 ; Vogel et al. 2009 ; Natali et al. in review ) Our results of higher NDVI in areas with higher roughness also support the idea of increased produc tivity as permafrost thaws (Figure 3 2B), as do several other studies that show NDVI is positively correlated with both R eco and GPP ( Boelman et al. 2003 ; Vourlitis et al. 2003 ; La Puma et al. 2007 ) Ecosystem respiration also increased in areas of the landscape with greater roughness d uring the GS of both years (Figure 3 5B & Figure 3 3) likely because of greater temperature and moisture a ssociated with greater ALD (Figure 3 2A; ( Vogel et al. 2009 ; Lee et al. 2010 ) ). More organic C is exposed to above freezing temperatures as ALD increases. These abiotic changes stimulate decomposition and nitrogen mineralization, which result in increased heterotrophic respiration ( Shaver et al. 1992 ) Vogel et al (2009) found that as subsidence incre ased, ALD increased and, in conjunction, both GPP and R eco increased. These results are
57 consistent with a large body of work showing that temperature and moisture are often the major determinates of organic matter decomposition and ecosystem respiration ( Oberbauer et al. 1991 ; Shaver et al. 1992 ; Hobbie et al. 2000 ; Xu & Qi 2001 ; Davidson et al. 2002 ; Davidson & Janssens 2006 ) There is also increased autotrophic respiratio n from more highly productive plants in subsided areas contributing to the overall increase in R eco during the GS ( Schuur et al. 2007 ; Vogel et al. 2009 ) In contrast to the GS, areas with increased roughness had lower C e mi ssions during the winter (Figure 3 5B & Figure 3 3). Even though areas with less roughness had higher R eco during winter, the overall C emission throughout 2008 and 2009 was still greater from areas with highest roughness. The reversal of the relationship of roughness and R eco during the winter is opposite previous work at the site that estimated more subsided areas have greater C emissions ( Vogel et al. 2009 ) They attributed greater winter C flux from subsided areas to warmer soils resulting from delayed active layer refreezing and the added insulation from snow accumulating in subsided areas, but data in the critical winter months was admittedly scarce ( Hinkel & Hurd 2006 ; Vogel et al. 2009 ) The inconsistency of our data may be due to differential diffusio n through variations in snow cover trapped in subsided areas. The coefficient of variation of the top GAM was also very low, only 0.16 and 0.22 during winter of 2008 and 2009, respectively. Overall, there was little variation in winter C flux (throughout s pace or time) and we believe more measurements are needed before our winter pattern can be fully supported. More winter change. Concluding Remarks We estimate d the C balance of a heterogeneous landscape undergoing permafrost thaw by incorporating spatial variation into an eddy covariance estimate. We found strong relationships
58 between thaw induced ground subsidence and ALD and NDVI, which both correlate with C flux. These microtopographical changes also strongly correlated with NEE. By using GAMs, we incorporated these spatial relationships back into final EC C balance estimates during gap filling. Thus, we achieved a more accurate C estimate for the heterogeneo us landscape and could make predictions for areas undergoing various degrees of permafrost thaw. Using GAMs, we were better able to predict C exchange during seasonal transitions, which indicates this type of gap filling strategy would be good in systems w ith high temporal variability. Because all natural ecosystems vary through space and time, we believe GAMs can be an important tool for achieving more accurate C estimates. The use of GAM's will also allow EC towers to be placed in more heterogeneous envir onments than they have been previously used. As permafrost thaws within this upland tundra ecosystem, a heterogeneous environment is created by changes in microtopography. We found this ecosystem was a C sink during two consecutive years, and areas with gr eater thaw exhibited greater C sequestration (GPP) and greater C loss (R eco ). Thawing of permafrost increases the amplitude of the C cycle, which has important implications for the future landscape level C balance ( Zimov et al. 1996 ) Currently, GPP is stimulated more than R eco but this balance may shift because we found that NDVI diminished with increased permafrost thaw indicating there may be an upper limit in productivity, unless successional changes in veget ation occur. Although we found the ecosystem was a C sink during the measurement campaigns of both years, this is not representative of the annual C balance. We did not measure C flux during a portion of the winter season, and even though C fluxes during t his time period are relatively low, the length of the season make it very important. By linearly extrapolating between these missing
59 winter periods, we found that annually the ecosystem became a C source of 60 gCm 2 year 1 and 13 gCm 2 year 1 in 2008 and 200 9, respectively.
60 Table 3 1 Coefficient of variation (R 2 mean square error (RMSE standard deviation) for GAM and NL models during GS day (GS D), GS nights (GS N), and winter of 2008 and 2009. Time P eriod Model R 2 RMSE GS D 2008 GAM 0.83 0 0.0670.03* NL 0.78 2458.6 0.1610.06 GS N 2008 GAM 0.26 0 0.1320.09 NL 0.26 93.4 0.1510.10 Winter 2008 GAM 0.16 0 0.0470.04 NL 0.01 116.9 0.0530.05 GS D 2009 GAM 0.83 0 0.0520.02* NL 0.78 5971.8 0.1810.17 GS N 2009 GAM 0.29 0 0.0850.06 NL 0.22 772.9 0.1360.09 Winter 2009 GAM 0.22 0 0.0410.03 NL 0.12 363.7 0.0650.06 Significantly different than model counterpart at p<0.05 Table 3 2 Carbon estimates (GPP, R eco NEE) in g C m 2 for the wind direction with the minimum roughness (Min), maximum roughness (Max) and the average (Avg) of all 360 wind directions during June to December 2008 and April to October 2009. Negative numbers denote when the ecosystem is taking up carbon. 20 08 2009 Roughness GPP R eco NEE GPP R eco NEE Min 300.3 288.3 11.2 403.6 386.6 33.6 Max 397.4 291.9 106.1 586.4 450.6 149.8 Avg 337.1 289.5 47.5 498.7 410.3 87.8
61 Figure 3 1. Map of microtopography and 360 transects surrounding the EC tower A) Map of microtopography surrounding the EC tower ( ), with lighter shades indicating areas where the ground surface is higher than the mean elevation of the landscape, and darker shades where the ground surface was lower than the mean elevation (i.e. subsided). B) 360 transects radiating out from EC tower ( ) corresponding to the wind direction sample by the tower. The color of transects grades from white to black as the degree of roughness increases. Note that in general the roughest transects o ccur to the North and North West of the EC tower, while transects to the South and East have lower roughness.
62 Figure 3 2 R elationship s betw een microtopography and ecosystem properties A) Non linear relationship between active layer depth (ALD; cm) and microtopography (adj. R 2 =0.51) with 95% confidence intervals. Note the small amount of variation in ALD where microtopography is positive or slightly negative, then an exponential increase in ALD as microtopography falls below 0.2 m. B) The non linear relationship between mean NDVI and roughness (adj. R 2 =0.55) with 95% confidence intervals. NDVI linearly increases until roughness reaches 0.14, then NDVI levels out or decreases.
63 Figure 3 3. Relationships between roughness and net ecosystem exchange (gC m 2 ) during GS nights, GS day and winter. Actual data points displayed behind the trend lines.
64 Figure 3 4. Carbon fluxes for the three methods of prediction. A) Gross primary production (GPP: gCm 2 ), B) ecosystem respiration (R eco : gCm 2 ), and C) net ecosystem exchange (NEE: gCm 2 ) for the three methods of prediction: GAM 360, GAM 1 and NL over the measurement campaigns of 2008 and 2009. Shaded areas denote winter conditions an d negative values occur when the ecosystem is a C sink.
65 Figure 3 5. Carbon fluxes for different wind direction s A) Gross primary production (GPP: gCm 2 ), B) ecosystem respiration (R eco : gCm 2 ), and C) net ecosystem exchange (NEE: gCm 2 ) for the landscape on average (360 Average) and for the wind direction with max roughness and min roughness over the measurement campaigns of 2008 and 2009. Shaded areas denote winter conditions and negative values occur when the ecosystem is a sink.
66 C HAPTER 4 TUNDRA ECOSYSTEMS OBSERVED TO BE CO2 SOURCES DUE TO DIFFERENTIAL AMPLIFICATION OF THE CARBON CYCLE Background Over the past five decades, researchers have tried to determine whether high latitude ecosystems will become a source of carbon (C) to t he atmosphere, or remain a sink, as climate changes. This question is of global importance because these ecosystems cover large areas of the northern high latitudes and store vast quantities of organic C in their soils ( Schuur et al. 2008 ; Tarnocai et al. 2009 ) Based on their large standing stocks of carbon, researchers have determined that these ecosystems have been, on average, a C sink for the past 10,000 yea rs ( Harden et al. 1992 ; Hicks Pries et al. 2 012 ) The ability of tundra to sequester and store carbon is due to long, harsh winters and poorly drained permafrost soils, which create conditions that slow decomposition relative to plant production ( Chapin et al. 1980 ; Miller et al. 1983 ; Billings 1987 ; Post 1990 ; Oechel & Billings 1992 ; Hobbie et al. 2000 ) But temperatures in high latitudes are rising ( Chapin et al. 2005 ; IPCC 2007 ) The response of the carbon cycl e to this climate forcing is of vital importance but the magnitude and timeline of change is uncertain. Two opposing feedbacks within the carbon cycle will determine the future C balance of tundra ecosystems. Rising temperatures resulting from increased at mospheric CO 2 concentrations could warm and thaw permafrost soils and stimulate decomposition and ecosystem respiration, resulting in a positive feedback to climate change by further increasing atmospheric CO 2 concentration (Shaver et al. 1992; Hobbie et a l. 2002; Schuur et al. 2008; McGuire et al. 2009; Grosse et al. 2011). On the other hand, increases in temperature and CO 2 concentration could stimulate primary production. Stimulation effects could be direct through effects on plant Belshe EF, EAG Schuur, BM Bolker (submitted) Tundra ecosystems observed to be CO 2 sources d ue to differential amplification of the carbon cycle. Ecology Letters.
67 physiology, or indirec t through increases in plant available nutrients released from decomposing organic matter and/or lengthening of the growing season ( Shaver et al. 1992 ; Johnson et al. 2000 ) If gross primary production e xceeds ecosystem respiration, the ecosystem will sequester CO 2 from the atmosphere and act as a negative feedback to climate change. Much effort has gone into quantifying the carbon balance of tundra ecosystems and understanding the controls over carbon up take and emission. Previous studies of growing season CO 2 exchange have shown initial release of CO 2 in the 1980s ( Oechel et al. 1993 ) followed by a longer term response of increased growing season CO 2 uptake ( Oechel et al. 2000 ; Ueyama et al. in press ) However, growing season patterns of CO 2 flux are only a part of the picture and by themselves do not provide complete information about the trajectory of tundra ecosystems. Winter CO 2 emissions have been recognized as an integral part of tundra carbon balance, and although there is considerably less winter data in the literature, recent efforts have improved our knowledge of the controls over winter CO 2 flux ( Fahnestock et al. 1998 ; Fahnestock et al. 1999 ; Welker et al. 2000 ; Grogan & Jonasson 2006 ; Nobrega & Grogan 2007 ; Sullivan et al. 2008 ; Rogers et al. 2011 ) However, a regional synthesis of winter CO 2 emission is missing. In this meta analysis we compiled CO 2 flux observations from both growing season and winter time periods. Using time series analysis, we assess whether seasonal or annual CO 2 fluxes have changed over time. By comparing CO 2 fluxes from different sites, we investigated the response of the carbon cycle to local differences in climate. Our intent is to re visit the question of whether tundra ecosystems are sources or sinks of CO 2 and to determine if we have observational evidence to support either trajectory.
68 Material and Methods Literature Search We compiled published observational data on CO 2 flux (Gross Primary Production: GPP, Ecosystem Respiration: ER, Net Ecosystem Exchange: NEE) from fifty four studie s at 32 different sites (Figure 4 1) ( Johnson 1970 ; Coyne 1975 ; Poole 1982 ; Grulke et al. 1990 ; Giblin et al. 1991 ; Oechel et al. 1993 ; Oechel et al. 1995 ; Oberbauer et al. 1996 ; Oechel et al. 1997 ; Vourlitis & Oechel 1997 ; Fahnestock et al. 1998 ; Hob bie 1998 ; Oberbauer et al. 1998 ; Oechel et al. 1998 ; Fahnestock et al. 1999 ; Grogan & Chapin 1999 ; Jones et al. 1999 ; Vourlitis & Oechel 1999 ; Welker et al. 1999 ; Oechel et al. 2000 ; Soegaard et al. 2000 ; Vourlitis et al. 2000a ; Vourlitis et al. 2000b ; Welker et al. 2000 ; Zamolodchikov et al. 2000 ; Heikkinen et al. 2002 ; Harazono et al. 2003 ; Epstein et al. 2004 ; Heikkinen et al. 2004 ; Grogan & Jonasson 2006 ; Kwon et al. 2006 ; Groendahl et al. 2007 ; Kutzbach et al. 2007 ; Lafleur & Humphreys 2007 ; Nobrega & Grogan 2007 ; van der Molen et al. 2007 ; Arens et al. 2008 ; Nobrega & Grogan 2008 ; Sjgersten et al. 2008 ; Sullivan et al. 2008 ; Vogel et al. 2009 ; Huemmrich et al. 2010 ; Taggeson et al. 2010 ; Zona et al. 2010 ; Humphreys & Lafleur 2011 ; Olivas et al. 2011 ; Rocha & Shaver 2011 ; Rogers et al. 2011 ; Belshe et al. 2012 ; Euskirchen et al. 2012 ; Trucco et al. 2012 ; Ueyama et al. in press ; Natali et al. in review ) NEE is the net exchange of CO 2 between the atmosphere and the ecosystem ove r a time interval ( Baldocchi 2003 ) and is the balance of two strong, opposing processes during the growing seaso n: CO 2 uptake by primary producers (gross primary production, GPP) and respiration losses of CO 2 by both primary producers and heterotrophs (ecosystem respiration, ER). During winter, tundra vegetation is covered by snow and photosynthesis is negligible, t herefore NEE is effectively equivalent to ER. We used the convention that negative values indicate CO 2 uptake by the ecosystem ( Baldocchi 2003 )
69 We gathered CO 2 flux data by conducting an Internet survey of ISI Web of Science and 2 flux observations from tundra ecosystems. As a second line of investigation, we contacted all lead investigators identified during our initial literature search and inquired about any additional studies that were previously overlooked. As a second filter for our meta analysis, we only included studies if flux measurements were distributed throughout the season (growing season or winter) and spanned the range of conditions encountered at each site. As a third filter, we only this category includes a variety of plant communities across a moisture gradient. We specifically excluded studies from high latitude forests, as well as fens, bogs, and mires, because these wetlands have different controls over carbon balance. We were in terested in non manipulated fluxes but also report control or ambient CO 2 flux estimates from experimental manipulations. For comparison, all data are reported as total seasonal fluxes in g C m 2 for either growing season or wi nter time periods Because th ere was a wide range in the interpretation of what constitutes the growing season and winter, and to follow conventions comparable with previous meta analyses ( Oechel et al. 1993 ; Oechel et al. 2000 ) we standardized the length of the growing season (100 days) and winter (245 days) by di viding seasonal estimates by the length of the reported season to estimate average daily (g C m 2 d 1 ) flux. We then multiplied daily estimates by the standardized growing season and winter length. We chose season lengths based on the average length of the seasons reported in our meta analysis, which is why the length of the growing season (100 days) and winter (245 days) does not add up to a full year. During the remaining 20 days, at the season transitions, CO 2 fluxes are roughly balanced with GPP=ER (i.e net flux is approximately zero), so excluding these periods does not alter estimates of annual flux derived by adding growing
70 season and winter estimates. Data were gathered on climate corresponding to the year when CO 2 fluxes were measured (mean annual temperature: MAT, total annual precipitation: TAP) for ea ch of the 32 sites sampled (Figure 4 1) from the University of Delaware precipitation and air temperature data set on NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, at http://www.esrl.noaa.gov/psd/. These weather data are a gridded interpolation of land station data and have a spatial resolution of 0.5 x 0.5 degrees, which may reduce data variation of near by sites. If data were reported as an average of multiple sites, we report the average climate c onditions. Data Analysis To elucidate changes in the tundra carbon cycle over space and time, we explored both temporal trends in CO 2 fluxes and relationships among CO 2 fluxes and spatial differences in climate. First, we examined how seasonal CO 2 fluxes have changed over the last four decades, and estimated the historical annual CO 2 balance from our temporal estimates. Temporal trends in growing season and winter CO 2 fluxes were estimated using a linear mixed effect model with site and site x year interaction as random effects to account for variation in the slopes with respect to time due to local site differences. There was a small but significant correlation between year and long term average (1950 2000) temperature (r=0.21; p =0.003), and long t erm average precipitation (r=0.32; p =0.001) for sites included in the analysis. An autoregressive (AR1) correlation structure was added to the model to account for temporal autocorrelation within sites. Flux estimates were aggregated by site for each year and weighted by the number of observations within each site year combination. Simulations of the null hypothesis (no overall trend at the population level) suggested that this approach was anticonservative, with a type I error rate of approximately 0.1 for a nominal p value of 0.05; therefore, we used parametric bootstrapping to achieve an appropriately conservative estimate for the significance of the trend
71 ( Booth 1995 ) Because our winter data were insufficient to support such a complex analysis, we did not include the correlation structure in the winter model and only included site (not the site x year inter action) as a random effect. To calculate an annual estimate of tundra CO 2 balance, we fitted the growing season model to a restricted data set containing only sites sampled during winter (to avoid calculating an annual estimate from disparate sites). We th en summed the predicted values from the restricted growing season and winter models to calculate the annual estimate. The variance of the annual estimate was estimated as the sum of the variances of the winter and growing season models ( Lyons 1991 ) test of the differences in estimated slopes in order to determine s ignificant seasonal differences in the rates of change through time. Second, we used simple linear regressions to explore if the observed variation in CO 2 fluxes (NEE, GPP, ER) during the growing season and winter were related to spatial differences in cli mate (MAT & TAP). Because MAT and TAP were highly correlated (r=0.7; p <0.001) for sites in our analysis, and the observational record of temperature is much stronger, we present only results from MAT, recognizing that responses to temperature and precipi t ation are strongly confounded We also explored the use of growing season (June Aug.) temperature to explain growing season trends in CO 2 fluxes, but found MAT explained more overall variation (higher adjusted R 2 ). In addition, slopes estimated with growing season temperature did not differ from MAT estimated slopes for either GPP or ER, and growing season temperature was a non significant predictor of growing season NEE ( p =0.47); therefore, we use MAT in our analysis of both seasons. To determine if the growing season or the winter was responding more strongly to test. Similar to the temporal analysis above, we calculated an annual response to temperature by adding
72 estimat ed growing season and winter trends together. All models were fit using the base and nlme packages, version 3.1 98 ( Pinheiro et al. 2011 ) and graphics were produced usi ng ggplot, version 0.8.9 ( Wickham 2009 ) in R (R Core Development Team 2011). Resul ts Changes in CO 2 Flux Over Time Over the past four decades, tundra ecosystems measured during the growing season have been found to be both CO 2 sources and CO 2 sinks, but overall growing season NEE has been decreasing over time (i.e., net CO 2 uptake is in creasing: Figure 4 2). Seventy percent of the total growing season NEE observations were below zero (net CO 2 uptake). Across the full dataset, growing season NEE has decreased 3.81.4 g C m 2 yr 1 (standard F test p =0.007; bootstrap adjusted p =0.03). Due to the limited number of sites and sparse observations at individual sites during the 1980s, we repeated this linear analysis using data collected post 1990, when the majority of observations were collected. The decreasing trend of growing season NEE for t his period is 4.50.9 g C m 2 yr 1 (bootstrap adjusted p = 0.005), which amounts to an estimated cumulative increase of 90 g C m 2 of net growing season C uptake between 1990 and 2010. Before 1990, the limited number of sites and sparse observations do not allow for robust tests of trends through time; thus the best fit line extending through the whole dataset is not well constrained for the early decades of the dataset. The trend in winter CO 2 emissions over time (1.51.1 g C m 2 yr 1 ), across the entire w inter dataset (1979 2010) was much weaker and not statistically signifi cantly different from zero (Figure 4 2B ; p =0.18), when including a random intercept to control for among site differences. However, the winter CO 2 flux dataset overall has many fewer s ites, and fluxes at individual sites are more sparsely quantified in comparison to growing season studies. Early records in the winter dataset came fr om relatively cold sites only. The middle period (1995
73 2000) contains a wider range of sites, but come mo stly from sites sampled only during one year. The most recent period of the record (from 2004 to 2010) contains sites that span the range of environmental conditions (MAT, MAP), and includes multiple sites monitored over several seasons that have more cohe rent with site t rends. In the latest period alone, we find a significant increase in winter CO 2 emissions over time, with an estimated increase of 15.94.3 g C m 2 yr 1 ( p =0.002;Figure 4 2B ). This increase in recent CO 2 emissions could be due to a recent amplification of winter CO 2 emissions or to the stronger data record, which improved our ability to estimate temporal trends in winter CO 2 flux. However, the exact slope of the line may not be representative of longer te rm trends because it is fit to only seven years of data. Estimates of annual CO 2 flux based on combining growing season and winter regressions show that the mean flux across the dataset ranged from an estimated net loss of 82 g C m 2 yr 1 in 1979 to an est imated loss of 21 g C m 2 yr 1 in 2010 (Figure 4 3). In the early 1980s it was not possible to determine if this net loss was statistically different from C neutral due to data scarcity. But for several decades from the mid 1980s until the early 2000s tund ra ecosystems were a net source of CO 2 to the atmosphere, with a mean net loss of 76 g C m 2 yr 1 in 1982 and a net loss of 36 g C m 2 yr 1 in 2002. By the mid 2000s, the overall declining but statistically insignificant annual trend (est= 2.01.9; p =0.18) pushed the estimated annual CO 2 flux towards net C neutrality. Clearly, observed increases in growing season CO 2 uptake are shifting tundra sites toward becoming a CO 2 sink on an annual basis. But the magnitude of the more uncertain winter emissions also plays a role. Using the winter trend from the recent data intensive period (2004 2010;Figure 4 3B ) counterbalances the shift toward C neutrality and reverses the trend through time, resulting in a continued increase in CO 2 em issi ons on an annual basis (Fig ure 4 3B ). Although it is hard to know how trends documented over a 7 year period will change in the
74 future, this highlights the important role of winter CO 2 emissions in determining annual CO 2 exchange, and emphasizes the need for more sustained winter CO 2 flux data from a variety of locations. Relationships Between CO 2 Flux and Spatial Differences in Climate To identify potential factors driving temporal trends in CO 2 flux, we explored patterns among CO 2 fluxes and spatial differences in MAT. We recognize that MAT is highly correlated with precipitation at the sites included in our analysis, so we are using MAT as a proxy for the general climate influence on CO 2 flux. Both growing season CO 2 emissions (Figure 4 4A ) and CO 2 uptake (Figure 4 4B ) are greater at sites with higher temperatures, with an overall negative trend between MAT and growing season NEE. This corresponds to higher net growing season CO 2 uptak e at sites with higher MAT (Figure 4 4C ). The slope of the growing season GPP trend is significantl y greater than the slope of the growing season ER trend (est= 3.62.0; p =0.05), suggesting that growing season GPP responds more strongly than growing season ER to differences in MAT. In addition, the analysis of growing season CO 2 flux by MAT, in contras t to the analysis with time, does not feature positive growing season CO 2 fluxes where CO 2 losses from respiration would need to be greater than plant uptake. This suggests that the period of net CO 2 loss predicted by the re gression model before 1990 (Figu re 4 2) was actually driven by data scarcity and influential points measured early in the record. During the winter, when only heterotrophic respiration is occurring, sites with higher MAT exhibited greater CO 2 emissions (Figure 4 4D ). The estimate of predicted net annual (winter plus growing season) CO 2 flux is positive throughout t he range of temperatures (Figure 4 5), yet there was no statistical difference between the absolute value of slopes of the growing season and winter trends (est.=1.92.0; p = 0.2). This suggests that both growing season and winter CO 2 flux have similar responses to changes in MAT, but the magnitude of the CO 2
75 emission exceeds the temperature response of CO 2 uptake. Based on predictions from the annual response, mean annual C em issions should range from 23 to 56 g C m 2 across the gradient of tundra temperatures, and a 1C increase in temperature would increase annual emissions by 2 g C m 2 These results indicate that the recent amplification of arctic temperature should have in creased both CO 2 uptake and CO 2 emission at a similar rate, but that ultimately the source strength of tundra ecosystems is greater on an annual basis. Discussion In spite of the spatial variation in climate, soil characteristics, vegetation composition, a nd site histories within the tundra biome ( Callaghan et al. 2004 ) our meta analysis detected an apparent regional amplification of the C cycle in recent decades. Growing season net CO 2 uptake has definitively increased since the 1990s, and trends point towards a potential increase in winter CO 2 emissions, esp ecia lly in t he last decade (Figure 4 3). Changing temperature, and the factors that covary with it, may have played a role in this amplification and based on the relationship with si te mean annual temperature (Figure 4 5), CO 2 emissions exceed CO 2 uptake across the ra nge of temperatures that occur in the tundra biome. This indicates a shift from the historic state of tundra as a CO 2 sink. Early analysis of Alaskan tundra reported a shift in growing season CO 2 balance from a historic C sink to a C source in response to warming and drying ( Oechel et al. 1993 ) As more flux data were compiled from Alaskan tundra, across site C O 2 balance was measured as net CO 2 uptake during the growing season. It was hypothesized that tundra ecosystems could have initially lost CO 2 in response to warming, then metabolically adjusted to continued warming while subsequently increasing CO 2 sink ac tivity ( Oechel et al. 2000 ) Our analysis using linear models does not allow us to quantify short term fluctuations in source and sink activity, and data limitations preclude us from quantifying growing season trends in the early part of the record.
76 While we do not r efute early data from Alaskan sites, we do not believe our pre 1990 growing season trend of CO 2 emissions adds any information useful for assessing this period. However, our larger post 1990 dataset clearly shows an increasing trend of growing season CO 2 u ptake (Figure 4 2A ). Our meta analysis over a longer time period and from 32 sites across the tundra biome shows a more gradual and continuing increase in CO 2 uptake during the growing season, on average, through the modern period. Similar to our finding o f increased net growing season CO 2 uptake, recent studies have reported that tundra was a growing season CO 2 sink during the past two decades ( McGuire et al. 2012 ; Ueyama et al. in press ) with increasing CO 2 uptake in a majority of tundra sites in the 2000s ( McGuire et al. 2012 ) In addition, long term (since the early 1980s) remote sensing observations show a greening Arctic ( Myneni et al. 1997 ; Jia & Epstein 2003 ; Nemani et al. 2003 ; Goetz et al. 2005 ; Sitch et al. 2007 ; Jia et al. 2009 ) and increased shrub encroachment ( Stow et al. 2004 ; Sturm et al. 2005 ) both suggesting a regional response of accelerated C uptake by tundra vegetation. By extrapolating our post 1990 growing season trend to the areal extent of tundra (10.5 x 10 6 km 2 ; ( McGuire et al. 1997 ; Callaghan et al. 2004 ) w e estimate that tundra on average sequestered 13780 Tg C during the growing seasons of the last few decades (1990 2006). Winter CO 2 trends were less clear in our analysis across the entire data set, but the last 7 years of the record show a strong increas e in CO 2 emissions (Fig ure 4 2B ). We acknowledge that the estimated slope from a 7 year trend may not match a decadal scale trend; based on the entire record it appears that the 7 year slope is likely an overestimate of longer trends (Fig ure 4 2B ). This va riability highlights the great need for additional sustained winter data to strengthen this record. Summing the growing season and longer term winter trend, we estimated that tundra
77 sites were CO 2 sources during the 1990s, but we were unable to differentia te tundra from carbo n neutral during the 2000s (Figure 4 3). Using only the more recent 7 year winter trend predicts increasing CO 2 emission annually for t he most recent time period (Figure 4 3B ), and the increasing winter and annual temperature trends (Fi gure 4 5) indicate that recent increases of arctic air temperature should lead to an increase in annual CO 2 emissions. Diverse lines of evidence support an amplification of the tundra C cycle in recent decades, but large uncertainty in past estimates have impeded the ability to accurately determine the annual C balance. Top down atmospheric inversion studies show the Arcti c as a C sink ( 410400 Tg C y 1 ) during the 1990s ( Baker e t al. 2006 ; McGuire et al. 2009 ) although recent inversion analyses were unable to distin guish Arctic tundra from C neutral in the 1990s ( 13; range= 321 to 140 Tg C yr 1 ) or the 2000s ( 117; range= 439 to 243 Tg C yr 1 ; ( McGuire et al. 2012 ) Retrospective analyses from process based models show an amplification of tundra C uptake and C emissions, but in the last few decades predict that tundra is a C sink. However, model predictions vary in sink strength and exhibit large temporal and spatial variability in source and sink activity ( Clein et al. 2000 ; McGuire et al. 2000 ; Sitch et al. 2003 ; Sitch et al. 2007 ; Grant et al. 2011 ) In a recent comparative analysis McGuire et al. (2012) found that summed estimates of annual C fluxes over the past few decades (1990 2006) from inversion models ( 96; range= 331 to 173 Tg C yr 1 ) were C neutral, whil e regional ( 177; range= 284 to 41 Tg C yr 1 ) and global ( 86; range= 205 to 1 Tg C yr 1 ) process based models predicted that tundra was a C sink, although sink strength could not be determined due to uncertainty in estimates ( McGuire et al. 2012 ) In addition, they reported that observationally based estimates of annual CO 2 balance over the last two decades ranged from C neutral (10; range= 10 to 28 Tg C yr 1 ) to a weak C sink ( 82; range= 134 to 30 Tg C yr 1 ; ( McGuire et al. 2012 ) In contrast
78 our mean annual estimate (Figure 4 3A ) extrapolated to the areal extent of tundra from the same time period predic ts tundra was a CO 2 source (462378 Tg C yr 1 ). The divergence of these observationally based estimates could be due to the larger variation in vegetation types and sampling duration and frequency of studies included in the McGuire et al (2012) meta analy sis. For tundra to be a historic C sink over centuries to millennia, C uptake had to be greater than C emissions across the range of conditions within the biome. Yet our annual temperature response based on the past several decades of measurements predicts that tundra sites are CO 2 sources across the range of temperatures (Fig ure 4 5). Within this range, tundra sites with warmer temperatures exhibit greater CO 2 fluxes (Figure 4 4; ( Ueyama et al. in press ) ) and both CO 2 uptake and CO 2 emissions are amplified by temperature a similar rate. These results suggest there has been a shift from the historic equilibrium and now C uptake and C emissions are of different magnitudes, with C emissions dominating across the range of conditions on a n annual basis. These predictions are supported by recent flux studies (2008 2010) from relatively warm (long term MAT of 1C) and cold (long term MAT of 12) tundra sites, which report CO 2 losses in five out of six annual observations ranging from 13 78 g C m 2 yr 1 and 2 82 g C m 2 yr 1 respectively ( Belshe et al. 2012 ; Euskirchen et al. 2012 ) A combina tion of several factors may have shifted the C balance of tundra ecosystems in recent decades. Surface air temperature has increased 0.35C per decade since the 1970s, with dramatic increases since the 1990s ( Polyakov et al. 2002 ; Hinzman et al. 2005 ; Serreze & Francis 2006 ; Euskirchen et al. 2007 ; McGuire et al. 2009 ) atmospheric CO 2 has increased by 50 ppm since 1980 ( Keeling & Whorf 2005 ) NOAA/ESLR), and the rate and areal extent of permafrost thaw has dramatically increased ( Jorgenson et al. 2001 ; Stow et al. 2004 ; Hinzman et al. 2005 ; Schuur et al. 2008 ; Osterkamp et al. 2009 ) This analysis indicates that tundra
79 ecosystems have not completely acclimated to recent changes and are currently undergoing a multi decade transition. In this transient state these ecosystems are CO 2 sources because of the unequal magnitude of CO 2 uptake and emissions. Based on our predicted annual temperature response, CO 2 emissions should exceed CO 2 uptake and push tundra to a CO 2 source on an annual basis, if they are not already doing so. Although our annual temporal trend is not definitive throughout the entire 40 year period, these data taken together suggest that despite increases in growing season CO 2 uptake, tundra ecosys tems are currently CO 2 sources on an annual basis.
80 Figure 4 1. Map of the World showing the 32 sites (stars) included in this meta analysis.
81 Figure 4 2. Temporal trends of net ecosystem CO 2 exchange (g C m 2 yr 1 ) (A ) During the growing season ( 3.81.4 g C m 2 yr 1 ; p =0.03). (B ) During the winter with the long term winter trend (1.51.1 g C m 2 yr 1 ; p =0.18) shown in solid blue and the most recent 7 year winter trend (15.94.3 g C m 2 yr 1 ; p =0.03) in dashed red. Estimates of slopes are report ed with standard errors; positive values denote a C source. The trend lines shown are from simple linear models, but slope estimates are based on the full mixed effect models.
82 Figure 4 3. Estimates of temporal trends of net ecosystem exchange (g C m 2 season 1 ) of CO 2 (A ) For the entire data record during the growing season (GS; blue), winter (red), and the annual period (black; produced by summing predicted values of seasonal trends and their varianc es). (B ) T emporal trends of GS, winter, and annual periods estimated from the last 7 years of the record. Trend lines are reported with 95% confidence intervals and are based on predictions from full mixed effect models; positive values denote a C source.
83 Figure 4 4. Linear relationships between mea n annual temperature and carbon fluxes. (A ) G rowing season ecosystem respiration (ER; est.=13.31.6; p =0.001; adj. R 2 =0.45). (B ) G rowing season gross primary production (GPP; est.= 16.72.3; p =0.001; adj. R 2 =0.43). (C ) G rowing season net ecosyst em exchange (NEE; est.= 4.51.7; p =0.007; adj. R 2 =0.04). (D ) W inter net ecosystem exchange (est.=6.81.3; p =0.001; adj. R 2 =0.31). Estimates of slopes are reported with standard errors; positive values denote a C source
84 Figure 4 5. Estimates of linear relationships between mean annual temperature (C) and net ecosystem exchange (g C m 2 season 1 ) of CO 2 during the growing season (GS; blue), and winter (red) with 95% confidence intervals, and the annu al estimate (black) of the temperature relationship produced by summing predicted values of seasonal trends and their variances. Positive values denote a C source.
85 CHAPTER 5 CONCLUSION Temperature is increasing in high latitude ecosystems and permafrost is dynamically shifting to be in equilibrium with contemporary climate ( Chapin et al., 2005, Jorgenson et al., 2006 Camill, 2005 Hinzman et al., 2005 Osterkamp et al., 2009 Romanovsky et al., 2010 Schuur et al., 2008 ) In response to warming soil active layers are thickening and thermokarst are forming, which creates a diversity of landforms that alter the physical environment in dynamic ways (Hinzman et al., 2005, Jorgenson and Osterkamp, 2005, Anisimov et al., 2007 Grosse et al., 2011). These physical changes feed back to a ffect biological and biochemical process es which ultimately alter the c ycling of carbon through high latitude ecosystems. Because of the large accumulation of soil carbon in this region and the potential for its release to the atmosphere ( Schuur et al., 2008 Tarnocai et al., 2009 ) it is imperative to develop ways to quantify both the spatial variation created by permafrost thaw and its effect on carbon fluxes. For this dissertation I focused on a particular type of thermokarst features (thermokarst gullies) that have formed within the Eight Mile Lake (EML) watershed and are typical and abundant in upland tundra landscapes. I detected and mapped thes e small, irregular thermokarst features by preforming a supervised cl assification on a high resolution IKONOS image I found that t welve percent of the EML watershed has undergone subsidence and thermokarst features predominantly form in valleys where tussock tundra resides. By 2008, 35% of the 3. 7 km 2 tussock tundra class had transitioned to thermokarst. These landscape level changes created by permafrost thaw at EML have important implications for its carbon cycle because thermokarst features are forming in soil carbon rich areas and are altering the hydrology in ways that increase seasonal thawing of the soil.
86 To better quantify how permafrost thaw effects the carbon balance of tundra at EML, I developed a modeling technique to incorporate spatial difference s in microtopography largely due to thermokarst subsidence into a carbon estimate. N et ecosystem exchange (NEE) of carbon was measured using eddy covariance (EC), during two 6 mo nth campaigns in 2008 and 2009, and a spatially explicit quantitative metric of permafrost thaw based on variation in microtopography was develo ped. This spatial variation was incorporated into an EC carbon flux estimate using a ge neralized additive model (GAM), which allowed us to make predictions about carbon exchange for the landscape as a whole, and for specific landscape patches throughout th e continuum of permafrost thaw and ground subsidence. In these study years, permafrost thaw appeared to increase the amplitude of the carbon cycle by stimulating both carbon release and sequestration, while the ecosystem remained a carbon sink at the lands cape scale. Although we found the ecosystem was a carbon sink during the measurement campaigns of both years, this is not representa tive of the annual carbon balance because we did not measure carbon flux during a portion of the winter season By linearly extrapolating between these missing winter periods, we found that annually the ecosystem became a carbon source of 60 gCm 2 year 1 and 13 gCm 2 year 1 in 2008 and 2009, respectively. To put these results into a broader context and determine if regionally tundra carbon balance has changed in recent decades I compile 40 years of CO 2 flux observations from 54 studies spanning 32 sites across northern high latitudes. Using time series analysis, I investigated if seasonal or annual CO 2 fluxes have changed over time, and whether spatial differences in mean annual temperature could help explain temporal changes in CO 2 flux. G rowing season net CO 2 uptake has definitely increased since the 1990s; the data also suggest (albeit less definitively) an increase in winter CO 2 emissions, especially in the last decade. In spite of the
87 uncertainty in the winter trend, tundra sites were annual CO 2 sources from the mid 1980s until the 2000s; and data from the last seven years show that tundra continue to emit CO 2 annually. In a ddition, CO 2 emissions exceed CO 2 uptake across the range of temperatures that occur in the tundra biome. Taken together, these data suggest that despite increases in growing season uptake, tundra ecosystems are currently CO 2 sources on an annual basis. Ov erall, the results from this work show that changes to permafrost are altering the landscape, and ultimately feeding back to amplify carbon c ycling within upland tundra in c entral Alaska. Although this work does not directly link permafrost degradation to the regional shift of tundra from a historic sink to a CO 2 source, the site level data points to this as a potential mechanism. Continued work is needed to fully understand and quantify the changes occurring to the carbon cycle as the landscape adjust to a warmer world. To better quantify annual carbon fluxes, more winter fluxes from a wider number of sites are needed This will allow for better annual estimates and enable the detect ion of temporal changes in winter fluxes. In addition to improved quantification of carbon fluxes, higher resolution digital elevation models and other spatial products are needed to facilitate the detection of the spatial variation created by permafrost thaw using remotely sensed images. Broadly this dissertation work tries to understand how climate change, and the spatial variation created by permafrost thaw, feed back to alter the carbon cycle in high latitude ecosystems. This is a difficult question because there are tremendous amount s of spa tial variation across the vast region underlain by permafrost in conjunction with spatial and temporal variation in changes in climate. My work has taken an empirical approach and by looking at patterns in observational data I have attempted to elucidate i f with current data we can detect changes occurring in the carbon cycle of tundra ecosystems. This work is in contrast to more mechanistic
88 studies that are focused on the determination of specific mechanisms that result in detected changes. It takes a comb ination of these two types of work to start to answer the difficult question of how the carbon cycle will change in a warmer world.
89 LIST OF REFERENCES Anderson D.R., Burnham K.P. & White G.C. (1998). Comparis on of Akaike information criterion and consistent Akaike information criterion for model selection and statistical inference from capture recapture studies. J Appl Stat 25, 263 282. Anderson D.R., Burnham K.P. & White G.C. (2001). Kullback Leibler informa tion in resolving natural resource conflicts when definitive data exist. Wildlife Soc B 29, 1260 1270. Anisimov O., Lobanov V., Reneva S., Shiklomanov N., Zhang T. & Nelson F. (2007). Uncertainties in gridded air temperature fields and effects on predicti ve active layer modeling. Journal of Geophysical Research Earth Surface 112, F02S14. Anisimov O. & Nelson F. (1996). Permafrost distribution in the Northern Hemisphere under scenarios of climatic change. Global and Planetary Change 14, 59 72. Arens S., S ullivan P.F. & Welker J.M. (2008). Nonlinear responses to nitrogen and strong interactions with nitrogen and phosphorus additions drastically alter the structure and Journal of Geophysical Research 113, G03S09. Aubinet M., Grelle A., Ibrom A., Rannik U., Moncrieff J., Foken T., Kowalski A.S., Martin P.H., Berbigier P., Bernhofer C., Clement R., Elbers J., Granier A., Grunwald T., Morgenstern K., Pilegaard K., Rebmann C., Snijders W., Valentini R. & Vesala T. (200 0). Estimates of the annual net carbon and water exchange of forests: The EUROFLUX methodology. Adv Ecol Res 30, 113 175. Aubinet M., Heinesch B. & Longdoz B. (2002). Estimation of the carbon sequestration by a heterogeneous forest: night flux corrections heterogeneity of the site and inter annual variability. Global Change Biology 8, 1053 1071. Baker D.F., Peylin P., Law R.M., Gurney K.R., Rayner P., Denning A.S., Bousquet P., Bruhwiler L., Chen Y.H., Ciasis P., Fung I., Heimann M., John J., Maki T., Ma ksyutov S., Masarie K., Prater M., Pak B., Taguchi S. & Zhu Z. (2006). TransCom 3 inversion intercomparison: impact of transport model errors on the interannual variability of regional CO2 fluxes, 1988 2003. Global Biogeochemical Cycles 20, 1002. Baldocch i D. (2003). Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future. Global Change Biology 9, 479 492. Belshe E.F., Schuur E.A.G., Bolker B.M. & Bracho R. (2012). Incorporating Spatial Heterogeneity Created by Permafrost Thaw into a Landscape Carbon Estimate. Journal of Geophysical Research: Biogeosciences 117. Beyer H.L. (2004). Hawth's Analysis Tools for ArcGIS. Billings W. (1987). Carbon balance of Alaskan tundra and tiaga ecosystems: past, present, and future. Quaternary Science Reviews 6, 165 177.
90 Boelman N., Stieglitz M., Rueth H., Sommerkorn M., Griffin K., Shaver G. & Gamon J. (2003). Response of NDVI, biom ass, and ecosystem gas exchange to long term warming and fertilization in wet sedge tundra. Oecologia 135, 414 421. Bolker B.M. (2010). bbmle: Tools for general maximum likelihood estimation. Bonan G.B. (2002). Ecological Climatology Cambridge University Press. Booth J. (1995). Bootstrap methods for generalized linear mixed models with applications to small area estimation. In: Statistical Modeling (ed. Seeber GUH, B.J. Francis, R. Hatzinger, G. Steckel Berger,). Springer New York, NY, pp. 43 51. Brown J. Ferrians O.J., Heginbottom J.A. & Melnikov E.S. (1998). International Permafrost Association circum Arctic map of permafrost and ground ice conditions. In. U.S. Geological Survey. Burba G.G., Anderson D.J., L. X. & McDermitt D.K. (2008). Addressing the influence of instrument surface heat exchange on the measurments of CO2 flux form open path gas analyzers. Global Change Biology 14, 1854 1876. Callaghan T., Bjorn L., Chernov Y., Chapin T., Christensen T., Huntley B., Ims R., Johansson M., Jolly D., Jona sson S., Matveyeva N., Panikov N., Oechel W. & Shaver G. (2004). Effects on the function of arctic ecosystems in the short and long term perspectives. Ambio 33, 448 458. Camill P. (2005). Permafrost thaw accelerates in Boreal peatlands during late 20th c entury climate warming. Climate Change 68, 135 152. Chapin F., Miller P.C., Billings W.D. & Coyne P.I. (1980). Carbon and nutrient budgets and their control in coastal tundra. In: An arctic tundra ecosystem. The coastal tundra at Barrow, Alaska. (eds. Bro wn J, Miller P, LL T & Bunnell F). Dowden Hutchinson and Ross Stroudsburg, Pennsylvania, pp. 458 482. Chapin F., Sturm M., Serreze M.C., McFadden J.P., Key J.R., Lloyd A.H., McGuire A.D., Rupp T.S., Lynch A.H., Schimel J.P., Beringer J., Chapman W.L., Epst ein H.E., Euskirchen E.S., Hinzman L.D., Jai G., Ping C.L., Tape K.D., Thompson C.D.C., Walker D.A. & Welker J.M. (2005). Role of Land Surface Changes in Arctic Summer. science 310, 657 660. Chapin F.S., Eugster W., McFadden J.P., Lynch A.H. & Walker D.A. (2000). Summer differences among arctic ecosystems in regional climate forcing. Journal of Climate 13, 12. Chapin F.S. & Shaver G. (1996). Physiological and growth responses of Arctic plants to a field experiment simulating climate change. Ecology 77, 8 22 840.
91 Clark K.L., Gholz H.L., Moncrieff J.B., Cropley F. & Loescher H.W. (1999). Environmental controls over net exchanges of carbon dioxide from contrasting Florida ecosystems. Ecological Applications 9, 936 948. Clein J.S., B. Kwiatkowski, A.D. McGuir e, J.E. Hobbie, E.B. Rastetter, J.M. Melillo & Kicklighter D.W. (2000). Modeling carbon responses of tundra ecosystems to historical and projected climate: a comparison of a plot and a global scale ecosystem model to identify process based uncertainties. Global Change Biology 6, 127 140. Congalton R.G. & Green K. (1999). Assessing the accuracy of remotely sensed data: Principles and practices Lewis Publishers, Boca Raton, Florida. Coyne P.I., J.J. Kelley (1975). CO2 exchange over the Alaskan arctic tundra: meteorological assessment by an aerodynamic method. Journal of Applied Ecology 12, 587 661. Davidson E., Belk E. & Boone R. (2002). Soil water content and temperature as independ ent or confounded factors controlling soil respiration in a temperate mixed hardwood forest. Global Change Biology 4, 217 227. Davidson E. & Janssens I. (2006). Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. Nature 440, 165 173. Davis N. (2001). Permafrost: A guide to frozen ground in transition. University of Alaska Press, Fairbanks, AK. Dioumaeva I., Trumbore S., Schuur E., Goulden M., Litvak M. & Hirsch A. (2003). Decomposition of peat from upland boreal forest: T emperature dependence and sources of respited carbon: Comparison of carbon exchange between boreal black spruce forest and the atmosphere for a wildfire age sequence. Journal of Geophysical Research 108, WFX3.1 WFX3.12. Dormann C.F. (2007). Effects of inc orporating spatial autocorrelation into the analysis of species distribution data. Global Ecology and Biogeography 16, 129 138. Epstein H.E., Beringer J., Gould W., Lloyd A.H., Thompson C.D., Chapin F.S., Michaelson G.L., Ping C.L., Rupp T.S. & Walker D.A (2004). The nature of spatial transitions in the Arctic. Journal of Biogeography 31, 1917 1933. Euskirchen E.S., Bret Harte M.S., Scott G.L., Edgar C. & Shaver G.R. (2012). Seasonal patterns of carbon dioxide and water fluxes in three representative tun dra ecosystems in northern Alaska. Ecosphere 3. Euskirchen S.E., McGuire A.D. & Chapin III F.S. (2007). Energy feedbacks of northern high latitude ecosystems to the climate system due to reduced snow cover during the 20th century warming. Global Change Bi ology 13, 2425 2438.
92 Fahnestock J., Jones M., Brooks P., Walker D.A. & Welker J.M. (1998). Winter and early spring CO2 efflux from tundra communities of northern Alaska. Journal of Geophysical Research 103, 29023 29027. Fahnestock J., Jones M. & Welker J .M. (1999). Wintertime CO2 efflux from Arctic soils: implications for annual carbon budgets. Global Biogeochemical Cycles 13, 775 779. Falge E., Baldocchi D., Olson R., Anthoni P., Aubinet M., Bernhofer C., Burba G., Ceulemans R., Clement R. & Dolman H. ( 2001). Gap filling strategies for defensible annual sums of net ecosystem exchange* 1. Agricultural and forest meteorology 107, 43 69. Giblin A., K. Nadelhoffer & G. Shaver (1991). Biogeochemical diversity along a riverside toposequence in arctic Alaska. Ecological Monographs 61, 415 435. Godin E. & Fortier D. (2012). Geomorphology of a thermo erosion gully, Bylot Island, Nunavut, Canada. Canadian Journal of Earth Science 49, 979 986. Goetz S., Bunn A.G., Fiske G.L. & Houghton R.A. (2005). Satellite obse rved photosynthetic trends across boreal North America associated with climate and fire disturbance. Proceedings of the National Academy of Sciences of the United States of America 102, 13521. Goulden M., Munger J., Fan S., Daube B. & Wofsy S. (1996). Mea surements of carbon sequestration by long term eddy covariance: methods and a critical evaluation of accuracy. Global Change Biology 2, 169 182. Goulden M.L., Wofsy S.C., Harden J., Trumbore S., Crill P.M., Gower S.T., Fries T., Daube B.C., S. M. F., Sutt on D.J., Bazzaz A. & Munger J.W. (1998). Sensitivity of boreal forest carbon balance to soil thaw. Science 279, 214 217. Goward S., Tucker C. & Dye D. (1985). North American vegetation patterns observed with the NOAA 7 advanced very high resolution radiom eter. Vegetation 64, 3 14. Grant R.F., Humphreys E.R., Lafleur P.M. & Dimitrov D.D. (2011). Ecological controls on net ecosystem productivity of a mesic arctic tundra under current and future climates. Journal Of Geophysical Research Biogeosciences 116, G01031. Groendahl L., Friborg T. & Soegaard H. (2007). Temperature and snow melt controls on interannual variability in carbon exchange in the high Arctic. Theoretical and Applied Climatology 88. Grogan P. & Chapin F.S. (1999). Arctic soil respiration: ef fects of climate and vegetation depend on season. Ecosystems 2, 451 459. Grogan P. & Jonasson S. (2006). Ecosystem CO2 production during winter in a Swedish subarctic region: the relative importance of climate and vegetation type. Global Change Biology 1 2, 1479 1495.
93 Grosse G., Harden J., Turetsky M., McGuire A.D., Camill P., Tarnocai C., Frolking S., Schuur E.A.G., Jorgenson T., Marchenko S., Romanovsky V., Wickland K.P., French N., Waldrop M., Bourgeau Chavez L. & Striegl R.G. (2011). Vulnerability of h igh latitude soil organic carbon in North America to distribance. Journal of Geophysical Research 116. Grulke N., G. Riechers, W. Oechel & U. Hjelm (1990). Carbon balance in tussock tundra under ambient and elevated atmospheric CO 2. Oecologia 83. Harazo no Y., M. Mano, A. Miyata & R. Zulueta (2003). Inter annual carbon dioxide uptake of a wet sedge tundra ecosystem in the Arctic. Tellus Series B Chemical And Physical Meteorology 55B, 215 231. Harden J., Mark R.K., Sundquist E.T. & Stallard R.F. (1992). D ynamics of soil carbon during deglaciation of the Laurentide ice sheet. Science 258, 1921 1924. Hastie T.J. & Tibshirani R.J. (1990). Generalized Additive Models Chapman and Hall/CRC, Boca Raton. Heikkinen J., V. Elsakov & Martikainem P.J. (2002). Carbon dioxide and methane dynamics and annual carbon balance in tundra wetland in NE Europe, Russia. Global Biogeochemical Cycles 16, 1115. Heikkinen J.E.P., T. Virtanen, J.T. Huttunen, V. Elsakov & Martikainen P.J. (2004). Carbon balance in East European tund ra. Global Biogeochemical Cycles 18, GB1023. Hicks Pries C.E., Schuur E.A.G. & Crummer K.G. (2012). Holocene carbon stocks and carbon accumulation rates altered in soils undergoing permafrost thaw. Ecosystems 15, 162 173. Hinkel K.M., Eisner W.R., Bockhe im J.G., Nelson F.E., Peterson K.M. & Dai X. (2003). Spatial extent, age, and carbon stocks in drained thaw lake basins on the Barrow Peninsula, Alaska. Arctic, Antarctic, and Alpine Research 35, 291 300. Hinkel K.M. & Hurd J.K. (2006). Permafrost destabilization and thermokarst following snow fence installation, Barrow, Alaska, USA. Arct Antarct Alp Res 38, 530 539. Hinzman L., Bettez N.D., Bolton W.R., Chapin F.S., Dyurgerov M.B., Fastie C.L., Griffith B., Hollister R.D., Hope A. & Huntington H.P. (2005). Evidence and implications of recent climate change in northern Alaska and other arctic regions. Climatic Change 72, 251 298. Hobbie S.E., F.S. Chapin III, (1998). The response of tundra plant biomass, aboveground pr oduction, nitrogen, and CO2 flux to experimental warming. Ecology 79, 1526 1544. Hobbie S.E., Schimel J.P., Trumbore S.E. & Randerson J.R. (2000). Controls over carbon storage and turnover in high latitude soils. Global Change Biology 6, 196 210.
94 Holling er D.Y., Kelliher F.M., Byers J.N., Hunt J.E., Mcseveny T.M. & Weir P.L. (1994). Carbon Dioxide Exchange between an Undisturbed Old Growth Temperate Forest and the Atmosphere. Ecology 75, 134 150. Huemmrich K.F., Kinoshita G., Gamon J.A., Houston S., Kwon H. & Oechel W.C. (2010). Tundra carbon balance under varying temperature and moisture regimes. Journal Of Geophysical Research Biogeosciences 115, G00I02. Humphreys E.R. & Lafleur P.M. (2011). Does earlier snowmelt lead to greater CO 2sequestration in tw o low Arctic tundra ecosystems? Geophysical Research Letters 38, L09703. IPCC (2007). Climate Change 2007: The Physical Science Basis. Contribution of Working Group 1 to the Fouth Assesment Report of the Intergovernmental Panel on Climate Change Cambridg e University Press, Cambridge, UK / New York, NY. Jafarov E.E., Marchenko S.S. & Romanovsky V.E. (2012). Numerical modeling of permafrost dynamics in Alaska using a high spatial resolution dataset. The Cryosphere Discussions 6, 89 124. Jensen J.R. (2005). Introductory digital image processing: A remote sensing perspective 3rd edn. Pearson Education, Inc., Upper Saddle River, NJ. Jia G. & Epstein H. (2003). Greening of arctic Alaska. 1981 2001. Geophysical Research Letters 30, 2067. Jia G., Epstein H. & Walker D.A. (2009). Vegetation greening in the canadian arctic related to decadal warming. Journal of Environmental Monitoring 11, 2231 2238. Johnson L., Shaver G., Cades D., Rastetter E., Nadelhoffer K., Giblin A., Laundre J. & Stanley A. (2000). Plant carbon nutrient interactions control CO2 exchange in Alaskan wet sedge tundra ecosystems. Ecology Johnson P.L., J.J. Kelley (1970). Dynamics of carbon dioxide and productivity in an arctic biosphere. Ecology 51, 73 80. Jones B.M., Grosse G., Arp C.D., Jo nes M.C., Walter Anthony K. & Romanovsky V. (2011). Modern thermokarst lake dynamics in the continuous permafrost zone, northern Seward Peninsula, Alaska. Journal of Geophysical Research 116, G00M03. Jones M., Fahnestock J. & Welker J.M. (1999). Early and late winter CO 2 efflux from arctic tundra in the Kuparuk River watershed, Alaska, USA. Arctic, Antarctic, and Alpine Research 31, 187 190. Jones M.C., Grosse G., Jones B.M. & Walter Anthony K. (2012). Peat accumulation in drained thermokarst lake basins in continuous, ice rich permafrost, northern Seward Peninsula, Alaska. Journal of Geophysical Research 117.
95 Jorgenson M. & Osterkamp T. (2005). Response of boreal ecosystems to varying modes of permafrost degradation. Canadian Journal of Forest Research 35, 2100 2111. Jorgenson M.T., Racine C.H., Walters J.C. & Osterkamp T.E. (2001). Permafrost Degradation and Ecological Changes Associated with a WarmingClimate in Central Alaska. Climatic Change 48, 551 579. Jorgenson M.T., Shur Y.L. & Pullman E.R. (2006). Abrupt increase in permafrost degradation in Arctic Alaska. Geophysical Research Letters 33. Kane D., Hinkel K., Goering D., Hinzman L. & Outcalt S. (2001). Non conductive heat transfer associated with frozen soils. Global and Planet ary Change 29, 275 292. Keeling C.D. & Whorf T.P. (2005). Atmospheric CO2 records from sites in the SIO air sampling network. In: In trends: A compendium of data on global change. Carbon dioxide information Analysis Center, Oak Ridge National Laboratory, U.S. Departmenf of Energy Oak Ridge, Tenn., U.S.A. Kim Y.J. & Gu C. (2004). Smoothing spline Gaussian regression: more scalable computation via efficient approximation. J Roy Stat Soc B 66, 337 356. Kljun N., Kormann R., Rotach N. & Meixner F. (2003). Com parison of the langrangian footprint model LPDM B with an analytical footprint model. Bound Lay Meteorol 106, 348 355. Kormann R. & Meixner F. (2001). An analytical footprint model for non neutral stratification. Boundary Layer Meteorol 99, 207 224. Kutz bach L., Wille C. & Pfeiffer E.M. (2007). The exchange of carbon dioxide between wet arctic tundra and the atmosphere at the Lena River Delta, Northern Siberia. Biogeosciences 4, 869 890. Kwon H. J., Oechel W.C., Zulueta R.C. & Hastings S.J. (2006). Effec ts of climate variability on carbon sequestration among adjacent wet sedge tundra and moist tussock tundra ecosystems. Journal Of Geophysical Research Biogeosciences 111, G03014. La Puma I., Philippi T. & Oberbauer S. (2007). Relating NDVI to ecosystem CO 2 exchange patterns in response to season length and soil warming manipulations in arctic Alaska. Remote Sensing of Environment 109, 225 236. Lafleur P.M. & Humphreys E.R. (2007). Spring warming and carbon dioxide exchange over low Arctic tundra in centra l Canada. Global Change Biology 14, 740 756. Laine A., Sottocornola M., Kiely G., Byrne K., Wilson D. & Tuittila E. (2006). Estimating net ecosystem exchange in a patterned ecosystem: example from blanket bog. Agricultural and Forest Meteorology 138, 231 243. Lantuit H., Pollard W.H., Couture N., Fritz M., Schirrmeister L., Meyer H. & Hubberten H .W. (2012). Modern and Late Holocene retrogressive thaw slump activity on the Yukon
96 Coastal Plain and Herschel Island, Yukon Territory, Canada. Permafrost and P eriglacial Processes 23, 39 51. Lantuit H. & Pollard W.H. (2008). Fifty years of coastal erosion and retrogressive thaw slump activity on Herschel Island, southern Beaufort Sea, Yukon Territory, Canada. Geomorphology 95, 84 102. Lawrence D.M., Slater A.G ., Romanovsky V.E. & Nicolsky D.J. (2008). Sensitivity of a model projection of near surface permafrost degradation to soil column depth and representation of soil organic matter. Journal of Geophysical Research 113. Lee H., Schuur E., Vogel J., Lavoie M. Bhadra D. & Staudhammer C. (2010). A spatially explicit analysis to extrapolate carbon fluxes in upland tundra where permafrost is thawing. Global Change Biology 17, 1379 1393. Liebethal C. & Foken T. (2007). Evaluation of six parameterization approaches for the estimation of ground heat flux. Theoretical and Applied Climatology 88, 43 56. Liebethal C., Huwe B. & Foken T. (2005). Sensitivity analysis for two ground heat flux calculation appro aches. Agricultural and Forest Meteorology 132, 253 262. Lyons L. (1991). A practical guide to data analysis for physical science students Cambridge University Press, Cambridge, England. Mack M., Schuur E., Bret Harte M., Shaver G. & Chapin III F. (2004) Ecosystem carbon storage in arctic tundra reduced by long term nutrient fertilization. Nature McGuire A., Melillo J.M., Kicklighter D.W., Pan Y., Xiao X., Helfrich J., Morre III M.B., Vorosmarty C.J. & Schloss A.L. (1997). Equilibrium responses of globa l net primary productivity and carbon storage to doubled atmospheric carbon dioxide: Sensitivity to changes in vegetation nitrogen concentration. Global Biogeochemical Cycles 11, 173 189. McGuire A.D., Anderson L.G., Christensen T.R., Dallimore S., Guo L., Haynes D.J., Heimann M., Lorenson T.D., MacDonald R.W. & Roulet N. (2009). Sensitivity of the carbon cycle in the Arctic to climate change. Ecological Monographs 79, 523 555. McGuire A.D ., Christensen T.R., Hayes D., Heroult A., Euskirchen E., Yi Y., Kimball J.S., Koven C., Lafleur P., Miller P.A., Oechel W., Peylin P. & Williams M. (2012). An assessment of the carbon balance of arctic tundra: comparisons among observations, proscess mode ls, and atmospheric inversions. Biogeosciences 9. McGuire A.D., Clein J.S., Melillo J.M., Kicklighter D.W., Meier R.A., Vorosmarty C.J. & Serreze M.C. (2000). Modelling carbon responses of tundra ecosystems to historical and projected climate: sensitivity of pan Arctic carbon storage to temporal and spatial variation in climate. Global Change Biology 6, 141 159.
97 Michaelson G.L., Ping C.L. & Kimble J.M. (1996). Carbon storage and distribution in tundra soils of Arctic Alaska, USA. Arctic Alpine Research 2 8. Mikan C., Schimel J. & Doyle A. (2002). Temperature controls of microbial respiration in arctic tundra soils above and below freezing. Soil Biology and Biochemistry 34, 1785 1795. Miller P., Kendall R. & Oechel W.C. (1983). Simulating carbon accumulati on in northern ecosystems. Simulation 40, 119 131. Moncrieff J.B., Massheder J.M., deBruin H., Elbers J., Friborg T., Heusinkveld B., Kabat P., Scott S., Soegaard H. & Verhoef A. (1997). A system to measure surface fluxes of momentum, sensible heat, water vapour and carbon dioxide. J Hydrol 189, 589 611. Morgenstern A., Grosse G., Gunther F., Fedorova I. & Schirrmeister L. (2011). Spatial analyses of thermokarst lakes and basins in Yedoma landscapes of the Lena Delta. The Cryosphere 5, 849 867. Myneni R. B., Keeling C.D., Tucker C.J., Asrar G. & Nemani R.R. (1997). Increased plant growth in the northern high latitudes from 1981 1991. Nature 386. Natali S., Schuur E. & Rubin R. (in review). Increased plant productivity in Alaskan tundra with experimental warming of soil and permafrost. Ecology Nemani R.R., Keeling C.D., Hashimoto H., Molly W.M., Piper S.C., Tucker C.J., Myneni R.B. & Running S.W. (2003). Climate driven increases in global terrestrial net primary production from 1982 1999. Science 300. Nobrega S. & Grogan P. (2007). Deeper snow enhances winter respiration from both plant associated and bulk soil carbon pools in birch hummock tundra. Ecosystems 10, 419 431. Nobrega S. & Grogan P. (2008). Landscape and ecosystem level controls on net carbon dioxide exchange along a natural moisture gradient in Canadian low arctic tundra. Ecosystems 11, 377 396. Oberbauer S., C. Gillespie, W. Cheng, A. Sala, R. Gebauer & J.D. Tenhunen (1996). Diurnal and seasonal patterns of ecosystem CO2 efflux from upland tundra in the foothills of the Brooks Range, Alaska, USA. Arctic and Alpine Research 28, 328 338. Oberbauer S., G. Starr & E.W. Pop (1998). Effects of extended growing season and soil warming on carbon dioxide and methane exchange of tussock tundra in Alaska. Journal of Geophysical Research 103, 29075 29082. Oberbauer S., Tenhunen J. & Reynolds J. (1991). Environmental effects on CO2 efflux from water track and tussock tundra in arctic Alaska, U.S.A. Arctic and Alpine Research 23, 162 169. Oechel W. & Billings W.D. (1992). Effects of global change on the carbon balance of arctic plants and ecosystems. In: Arctic ecosystems in a changing climate: an ecophysiological
98 per spective. (eds. Chapin III F, Jefferies R, Reynolds J, Shaver G & Svoboda J). Academic Press San Diego, CA, pp. 139 168. Oechel W., G. Vourlitis & S.J. Hastings (1997). Cold season CO2 emission from arctic soils. Global Biogeochemical Cycles 11, 163 172. Oechel W., Hastings S.J., Vourlitis G., Jenkins M., Riechers G. & Grulke N. (1993). Recent change of Arctic tundra ecosystems from a net carbon dioxide sink to a source. Nature 361, 520 523. Oechel W., Vourlitis G. & Hastings S. (1995). Change in Arctic C O^ 2 Flux Over Two Decades: Effects of Climate Change at Barrow, Alaska. Ecological Applications 5, 846 855. Oechel W., Vourlitis G., Hastings S., Ault R.P. & Bryant P. (1998). The effects of water table manipulation and elevated temperature on the net CO 2 flux of wet sedge tundra ecosystems. Global Change Biology 4, 77 90. Oechel W.C., Vourlitis G.L., Hastings S.J., Zulueta R.C., Hinzman L. & Kane D. (2000). Acclimation of ecosystem CO2 exchange in the Alaskan Arctic in response to decadal climate warmin g. Nature 406, 978 981. Olivas P., Oberbauer S., Tweedie C., Oechel W.C., Lin D. & Kuchy A. (2011). Effects of Fine Scale Topography on CO 2 Flux Components of Alaskan Coastal Plain Tundra: Response to Contrasting Growing Seasons. Arctic, Antarctic, and A lpine Research 43, 256 266. Osterkamp T. (2005). The recent warming of permafrost in Alaska. Global and Planetary Change 49, 187 202. Osterkamp T.E., Jorgenson M.T., Schuur E.A.G., Shur Y.L., Kanevskiy M.Z., Vogel J.G. & Tumskoy V.E. (2009). Physical and Ecological Changes Associated with Warming Permafrost and Thermokarst in Interior Alaska. Permafrost and Periglacial Processes 20, 235 256. Pinheiro J., Bates D., Sarkar D. & R development Core Team (2011). nlme: Linear and Nonlinear Mixed Effects Models R package version 3.1 98 edn. Polyakov I.D., Alekseev G.V., Bekryacy R.V., Bhatt U., Colony R., Johnson M.A., Karklin V.P., Makshtas A.P., Walsh D. & Yulin Y.V. (2002). Observationally based assessment of polar amplifica tion of global warming. Geophysical Research Letters 29, 1878. Poole D., P.C. Miller, (1982). Carbon dioxide flux from three Arctic tundra types in north central Alaska, USA. Arctic and Alpine Research 14, 27 32. Post W. (1990). Report of a workshop on c limate feedbacks and the role of peatlands, tundra, and boreal ecosystems in the global carbon cycle. In. Oak Ridge National Lab Oak Ridge, TN.
99 R Development Core Team (2010). R: A language and enviroment for statistical computing. URL http://www.R project.org/ Raynolds M.K., Walker D.A. & Maier H.A. (2006). NDVI patterns and phytomass distribution in the circumpolar Arctic. Remote Sensing of Environment 102, 271 281. Reynolds O. (1985). On the dy namical theory of incompressible viscous fluids and the determination of the criterion. Philosophical Transactions of the Royal Society A 186, 123 164. Richardson A., Mahecha M., Falge E., Kattge J., Moffat A., Papale D., Reichstein M., Stauch V., Braswel l B. & Churkina G. (2008). Statistical properties of random CO2 flux measurement uncertainty inferred from model residuals. agricultural and forest meteorology 148, 38 50. Rivkina E., Friedmann E., McKay C. & Gilichinsky D. (2000). Metabolic activity of p ermafrost bacteria below the freezing point. Applied and Environmental Microbiology 3239 3233. Rocha A. & Shaver G.R. (2011). Burn severity influences postfire CO2 exchange in arctic tundra. Ecological Applications 21, 477 489. Rogers M.C., Sullivan P.F. & Welker J.M. (2011). Evidence of Nonlinearity in the Response of Net Ecosystem CO 2 Exchange to Increasing Levels of Winter Snow Depth in the High Arctic of Northwest Greenland. Arctic, Antarctic, and Alpine Research 43, 95 106. Romanovsky V., Gruber S. Instanes A., Jin H., Marchenko S., Smith S., Trombotto D. & Walter K. (2007). Frozen Ground. In: Global outlook for ice and snow (ed. UNEP/GRID) Arendal, Norway. Romanovsky V.E., Smith S.L. & Christiansen H.H. (2010). Permafrost thermal state in the pola r Northern Hemisphere during the international polar year 2007 2009: a synthesis. Permafrost and Periglacial Processes 21, 106 116. Saito K., Kimoto M, Zhang T, Takata K & Emori S (2007). Evaluating a high resolution climate model: Simulated hydrothermal regimes in frozen ground regions and their change under the global warming scenario. J. Geophys. Res 112. Schmid H. (1997). Experimental design for flux measurements: matching scales of observations and fluxes. Agricultural and Forest Meteorology 87, 179 200. Schmid H. (2002). Footprint modeling for vegetation atmosphere exchange studies: a review and perspective. Agricultural and Forest Meteorology 113, 159 183. Schmid H. & Lloyd C. (1998). Spatial representativeness and the location bias of flux footpr ints over inhomogeneous areas. Agricultural and Forest Meteorology 93, 195.
100 Schuur E., Bockheim J., Canadell J.G., Euskirchen E., Field C.B., Goryachkin S.V., Hagemann S., Kuhry P., Lafleur P.M. & Lee H. (2008). Vulnerability of permafrost carbon to climate change: Implications for the global carbon cycle. BioScience 58, 701 714. Schuur E., J Vogel, K Crummer, H Lee, JO Sickman & TE Osterkamp (2009). The effect of permafrost thaw on old carbon rel ease and net carbon exchange from tundra. Nature 459, 556 559. Schuur E.A.G., Crummer K.G., Vogel J.G. & Mack M.C. (2007). Plant species composition and productivity following permafrost thaw and thermokarst in Alaskan tundra. Ecosystems 10, 280 292. Sel lers P. (1987). Canopy reflectance, photosynthesis and transpiration. II. The role of biophysics in the linearity of their interdependance. International Journal of Remote Sensing 6, 1335 1372. Serreze M. & Francis J.A. (2006). The arctic amplification de bate. Climate Change 76, 241 264. Shaver G., Billings W.D., Chapin III F.S., Giblin A.E., Nadelhoffer K.J., Oechel W.C. & Rastetter E.B. (1992). Global change and the carbon balance of arctic ecosystems. BioScience 42, 433 441. Shur Y., Hinkel K.M. & Nel son F.E. (2005). The transient layer: implications for geocryology and climate change science. Permafrost and Periglacial Processes 16, 5 17. Shur Y.L. & Jorgenson M.T. (2007). Patterns of permafrost formation and degradation in relation to climate and ec osystems. Permafrost and Periglacial Processes 18, 7 19. Sitch S., McGuire A.D., Kimball J., Gedney N., Gamon J., Engstrom R., Wolf A., Zhuang Q., Clein J. & McDonald K.C. (2007). Assessing the carbon balance of circumpolar Arctic tundra using remote sens ing and process modeling. Ecological Applications 17, 213 234. Sitch S., Smith B., IPrentice I.C., Arneth A., Bondeau A., Cramer W., Kaplans J.O., Levis S., Lucht W., Sykes M.T., Thonicke K. & Venevsky S. (2003). Evaluation of ecosystem dynamics, plant ge ography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biology 9, 161 185. Sjgersten S., van der Wal R. & Woodin S.J. (2008). Habitat type determines herbivory controls over CO2 fluxes in a warmer Arctic. Ecology 89, 2103 2116. Soegaard H., Nordstroem C., Frigorg T. & al. e. (2000). Trace gas exchange in a high arctic valley. Integrating and scaling CO2 fluxes from canopy to landscape using flux data, footprint modeling, and remote sensing. Global Biogeochemical Cycles 14, 725 744. Soil Survey Staff (1999). Soil Taxonomy: A Basic System of Soil Classification for Making and Interpreting Soil Surveys. In: Agricultural Handbook Washington, D.C.
101 Stephens J. L., Boggs K., G aribaldi A., Grunblatt J. & Helt T. ( 2001 ) Denali National Park and Preserve landcover mapping project Volume 1: Remote sensing data, procedures and results. In: National Resource Technical Report NPS/DENA/NRTR. Fort Collins, Colorado: National Park Service. Stieglitz M., Dery S.J., Romanovsky V.E. & Osterkamp T.E. (2003). The role of snow cover in the warming of Arctic permafrost. Geophysical Research Letters 30. Stow D., Hope A., McGuire D., Verbyla D., Gamon J., Huemmrich F., Houston S., Racine C., Sturm M., Tape K., Hinzman L., Yoshikawa K., Tweedie C., Noyle B., Silapaswan C., Douglas D., Griffith B., Jia G., Epstein H., Walker D., Daeschner S., Petersen A., Zh ou L. & Myneni R. (2004). Remote sensing of vegetation and land cover change in Arctic tundra ecosystems. Remote Sensing Environment 89, 281 308. Sturm M., Schimel J., Michaelson G., Welker J.M., Oberbauer S.F., Liston G.E., Fahnestock J. & Romanovsky V.E (2005). Winter biological processes could help convert arctic tundra to shrubland. BioScience 55, 17 26. Sullivan P., Welker J.M., Arens S.J.T. & Sveinbjornsson B. (2008). Continuous estimates of CO2 efflux from arctic and boreal soils during the snow c overed season in Alaska. Journal of Geophysical Research 113, G04009. Taggeson T., Mastepanov M., Tamstorf M.P., Eklundh L., Schubert P., Ekberg A., Sigsgaard C., Christensen T.R. & Strom L. (2010). Satellites reveal an increase in gross primary productio n in a greenlandic high arctic fen 1992 2008. Biogeosciences Discussions 7, 1101 1129. Tarnocai C., Canadell J.G., Schuur E.A.G., Kuhry P., Mazhitova G. & Zimov S. (2009). Soil organic carbon pools in the northern circumpolar permafrost region. Global Bio geochemical Cycles 23, 1 11. Thornley J.H. & Johnson I.R. (1990). Plant and Crop Modeling: A Mathmatical Approach to Plant and Crop Physiology Charldon, Oxford, U.K. Trucco C., Schuur E.A.G., Natali S., Belshe E.F., Bracho R. & Vogel J. (2012). Long term trends of CO2 exchange in a tundra ecoystem affected by permafrost thaw. Journal of Geophysical Research Biogeosciences 117. Tucker C. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing Environ 8, 127 150. Ueyama M., Iwata H., Harazono Y., Euskirchen E.S., Oechel W.C. & Zona D. (in press). Growing season and spatial variations of carbon fluxes of arctic and boreal ecosystems in Alaska. Ecological Applications van der Molen M.K., van Huissteden J., Parmentier F.J.W. & et al. (2007). The growing season greenhouse gas balance of a continental tundra site in the Indigirka lowlands, NE Siberia. Biogeosciences 4.
102 Venables W. & Ripley B. (2002). Modern Applied Statistics with S Fourth edn. Springer, New York, NY. Vogel J., Schuur E., Trucco C. & Lee H. (2009). Response of CO2 exchange in a tussock tundra ecosystem to permafrost thaw and thermokarst development. Journal of Geophysical Research 114. Vourlitis G. & Oechel W. (1997). Land scape scale CO 2, H 2 O vapour and energy flux of moist wet coastal tundra ecosystems over two growing seasons. Journal of Ecology 85, 575 590. Vourlitis G., Verfaillie J., Oechel W., Hope A., Stow D. & Engstrom R. (2003). Spatial variation in regional CO 2 exchange for the Kuparuk River Basin, Alaska over the summer growing season. Global Change Biology 9, 930 941. Vourlitis G.L., Harazono Y., Oechel W.C., Yoshimoto M. & Mano M. (2000a). Spatial and temporal variations in hectare scale net CO 2 flux, resp iration and gross primary production of Arctic tundra ecosystems. Functional Ecology 14, 203 214. Vourlitis G.L. & Oechel W.C. (1999). Eddy covariance measurements of CO2 and energy fluxes of an Alaskan tussock tundra ecosystem. Ecology 80, 686 701. Vour litis G.L., Oechel W.C., Hope A., Stow D., Boynton B., Verfaillie Jr J., Zulueta R. & Hastings S.J. (2000b). Physiological models for scaling plot measurements of CO2 flux across an arctic tundra landscape. Ecological Applications 10, 60 72. Wahrhaftig C. (1958). Quarternary and engineering geology in the central part of the Alaska Range. Geological Survey Professional Paper 293. Walker D. A., R aynolds M. K., D aniels F. J. A., Einarsson E., E lvebakk A., G ould W. A., K atenin A. E., Kholod S. S., M arkon C. J., M elnikov E. S., M oskalenko N. G., Talbot S. S., Y urtsev B. A. & other members of the CAVM team. ( 2005 ) The Circumpolar Arctic vegetation map. Journal of Vegetation Science, 16 267 282 Walter K.M., Smith S.L. & Chapin F.S. (2007). Methane bubbling f rom northern lakes: present and future contributions to the global methane budget. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 365, 1657 1676. Wang J., Sheng Y., Hinkel K.M. & Lyons L. (2012). Draine d thaw lake basin recovery on the western Arctic Coastal Plain of Alaska using high resolution digital elevation models and remote sensing imagery. Remote Sensing of Environment 119, 325 336. Webb E.K., Pearman G.I. & Leuning R. (1980). Correction of Flux Measurements for Density Effects Due to Heat and Water Vapor Transfer. Quarterly Journal of the Royal Meteorological Society 106, 85 100.
103 Welker J., Fahnestock J. & Jones M.H. (2000). Annual CO2 flux in dry and moist arctic tundra: field responses to inc reases in summer temperatures and winter snow depth. Climatic Change 44. Welker J.M., Brown K.B. & Fahnestock J. (1999). CO 2 flux in arctic and alpine dry tundra: comparative field responses under ambient and experimentally warmed conditions. Arctic, Ant arctic, and Alpine Research 31, 272 277. Wickham H. (2009). Springer, New York. Wilczak J., Oncley S. & Stage S. (2001). Sonic anemometer tilt correction algorithms. Boundary Layer Meteorol 99, 127 150. Williams M., Bell R., Spadavecchia L., Street E. & Van Wijk M. (2008). Upscaling leaf area index in an Arctic landscape through multiscale observations. Global Change Biology 14, 1517 1530. Wood S. (2006). Generalized Additive Models: An Introduction with R Chapman and Hall/CRC, Boca Raton. Wood S.N. (2008). Fast stable direct fitting and smoothness selection for generalized additive models. Journal of the Royal Statistical Scoiety (B) 70(3), 495 518. Xu M. & Qi Y. (2001). Soil surface CO2 efflux and its spatial and temporal variations in a young ponderosa pine plantation in northern California. Global Change Biology 7, 667 677. Yocum L., G.W. Adema & C.K. Hults (2006). A baseline study of permafrost in the T otlat Basin, Denali National Park and Preserve. In, pp. 37 40. Yoshikawa K. & Hinzman L.D. (2003). Shrinking thermokarst ponds and groundwater dynamics in discontinuous permafrost near Council, Alaska. Permafrost and Periglacial Processes 14. Zamolodchiko v D., Karelin D. & Ivaschenko A. (2000). Sensitivity of tundra carbon balance to ambient temperature. Water, Air, And Soil Pollution 119, 157 169. Zhang T., Barry RG, Knowles K, Heginbottom JA & Brown J (1999). Statistics and characteristics of permafrost and ground ice distribution in the Northern Hemisphere. Polar Geography 23, 132 154. Zhang T., Frauenfeld O., Serreze M., Etringer A., Oelke C., McCreight J., Barry R., Gilichinsky D., Yang D. & Ye H. (2005). Spatial and temporal variability in active la yer thickness over the Russian Arctic drainage basin. J. Geophys. Res 110. Zimov S., Davidov S., Voropaev Y., Prosiannikov S., Semiletov I., Chapin M. & Chapin F. (1996). Siberian CO2 efflux in winter as a CO2 source and cause of seasonality in atmospheri c CO2. Climatic Change 33, 111 120.
104 Zona D., Oechel W.C., Peterson K.M., Clements R.J., Paw K.T. & Ustin S.L. (2010). Characterization of the carbon fluxes of a vegetated drained lake basin chronosequence on the Alaskan Arctic Coastal Plain. Global Change Biology 16, 1870 1882. Zuur A., Ieno E., Walker N., Saveliev A. & Smith G. (2009). Mixed Effects Models and Extensions in Ecology with R Springer, New York, NY.
105 BIOGRAPHICAL SKETCH Fay Belshe grew up a poor Florida cracker and spent most of the days of her youth outdoors. Mostly, she was on the beach or in the shade of a large oak, trying to survive the summer heat. During this time she began to appreciate living things and wondered what created the beautiful patterns of na ture. Although she would be lying if she said that from that point on her life was dedicated to ecological sciences A lot of her years were spent drifting along, barely staying out of trouble. She finally became serious about school in the final years of her bachelor degree, when she studied horticultural sciences. After graduation, Fay continued studying plants in their natural environment and became interested in carbon cycling. During examined the spatial and temporal variation of seagr ass photosynthesis which included many fun summers snorkeling in the Florida Keys For her doctoral research, she expanded upon her knowledge of photosynthesis and focused on whole ecosystem carbon cycling. Specifically, understanding and quantifying land scape patterns in carbon cycling created as permafrost thaws in an upland tundra ecosystem. It would be a second lie if she said she enjoyed working in high latitude ecosystems, but she survived and got some pretty good data. Fay really enjoy s the intricat e patterns and stories that emerge from ecological data. By helping to tell the complex story of carbon, she hopes to contribute to the management of ecosystems in a warmer world. But beyond increasing knowledge within the scientific community, she hopes t o tell a story about our understanding of the earth system that can reach the general public, and can change their conception of the importance of carbon in our climate