The Polar Vegetation Photosynthesis and Respiration Model : a parsimonious , satellite-data-driven model of high-latitude CO 2 exchange

We introduce the Polar Vegetation Photosynthesis and Respiration Model (PolarVPRM), a remote-sensingbased approach for generating accurate, high-resolution (≥ 1 km, 3 hourly) estimates of net ecosystem CO2 exchange (NEE). PolarVPRM simulates NEE using polarspecific vegetation classes, and by representing high-latitude influences on NEE, such as the influence of soil temperature on subnivean respiration. We present a description, validation and error analysis (first-order Taylor expansion) of PolarVPRM, followed by an examination of per-pixel trends (2001–2012) in model output for the North American terrestrial region north of 55 N. PolarVPRM was validated against eddy covariance (EC) observations from nine North American sites, of which three were used in model calibration. Comparisons of EC NEE to NEE from three models indicated that PolarVPRM displayed similar or better statistical agreement with eddy covariance observations than existing models showed. Trend analysis (2001–2012) indicated that warming air temperatures and drought stress in forests increased growing season rates of respiration, and decreased rates of net carbon uptake by vegetation when air temperatures exceeded optimal temperatures for photosynthesis. Concurrent increases in growing season length at Arctic tundra sites allowed for increases in photosynthetic uptake over time by tundra vegetation. PolarVPRM estimated that the North American high-latitude region changed from a carbon source (2001–2004) to a carbon sink (2005–2010) to again a source (2011–2012) in response to changing environmental conditions.


Introduction
Large uncertainties presently exist in process-based model estimates of high-latitude North American NEE (Fisher et al., 2014), and limit understanding and monitoring of recent changes in the polar carbon cycle.Simultaneously, recent successes in generating site-level, data-driven estimates of net ecosystem CO 2 exchange (NEE) at Arctic sites (e.g., Shaver et al., 2007Shaver et al., , 2013;;Stoy et al., 2009), with little intersite variability in parameters (Loranty et al., 2011), have indicated the tremendous potential that exists for accurate estimates of regional-scale Arctic NEE to be modeled diagnostically from satellite observations.
In this article, we describe, validate and examine output from the newly developed Polar Vegetation Photosynthesis and Respiration Model (PolarVPRM).PolarVPRM is an arctic-specific remote-sensing-based model which estimates high-latitude NEE at a fine resolution (3 hourly, ≥ 1 km 2 ) using a diagnostic approach.PolarVPRM is presently in active use by the Arctic research community for a range of applications, including the examination and scaling-up of circumpolar eddy covariance (EC) observations of NEE, and as a priori estimates of Alaskan NEE for Lagrangian modeling of aircraft CO 2 concentration observations (Miller and Dinardo, 2012).

PolarVPRM formulation
PolarVPRM presents a high-latitude formulation of VPRM (Mahadevan et al., 2008).Both PolarVPRM and VPRM were Published by Copernicus Publications on behalf of the European Geosciences Union.
written in R (R Development Core Team, 2011), and provide straightforward yet effective calculations of terrestrial biospheric carbon exchange from remote-sensing observations.In both VPRM and PolarVPRM, NEE is calculated as the sum of respiration (R) and gross ecosystem exchange (GEE; the light-dependent portion of NEE), using the sign convention where CO 2 efflux to the atmosphere via R is positive, and CO 2 uptake through photosynthesis (GEE) is negative: (1) Relative to VPRM, PolarVPRM uses different inputs (described in Appendix A), vegetation classes (presented in Luus et al., 2013b) and model structure (described in Luus et al., 2013a), in order to ensure suitability for modeling high-latitude NEE.VPRM has previously been applied and validated across the USA and southern Canada (30-56 • N) (Mahadevan et al., 2008;Lin et al., 2011), and PolarVPRM is now applied to generate estimates of NEE across highlatitude regions (e.g., north of 55 • N).

Gross ecosystem exchange
GEE, or the photosynthetic uptake of C by vegetation, is calculated according to remote-sensing-based estimates of incoming shortwave radiation (SW; expressed as photosynthetically active radiation (PAR), where PAR = 1.98 × SW; Lin et al., 2011), air temperature (T air ), land surface water index (LSWI) from Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance, and estimates of the fraction of PAR absorbed by photosynthetically active vegetation (FAPAR PAV ), as estimated from the MODIS Enhanced Vegetation Index (EVI).
In Arctic tundra regions, GEE is therefore implicitly limited during the snow season, when EVI is decreased, suggesting that negligible amounts of photosynthetically active vegetation persist above the snow surface.Similarly, GEE can be limited when air temperatures are sub-optimal, or when vegetation is at an underdeveloped phenological stage.These limitations are implemented through use of dimensionless scaling variables T scale and P scale , respectively. (3) λ refers theoretically to the maximum light-use efficiency (LUE), or quantum yield, at low-light levels, but functions in practice as a combined LUE and scaling parameter.PAR 0 is the half-saturation value of PAR (Mahadevan et al., 2008;Lin et al., 2011).As in Mahadevan et al. (2008), T max = 40 • C and T min = 0 • C for all vegetation classes, and T opt = 20 • C for non-arctic vegetation classes.For barren/wetland regions (which include the Canadian High Arctic), a T opt = 10 • C, whereas T opt = 15 • C over shrub tundra and graminoid tundra, as approximated from literature (e.g., Tieszen, 1973;Chapin III, 1983).Plots of air temperature and growing season NEE at calibration sites were checked to ensure that these values appeared reasonable, but no optimization took place, to avoid correlation and instability of parameters (Mahadevan et al., 2008).
LSWI max refers to the maximum annual pixel-specific LSWI.LSWI is implemented as a limitation on GEE (W scale ) for forested regions north of 55 • N, just as in VPRM.However, water availability does not play a clear role in determining Arctic plant productivity (Oberbauer and Miller, 1979;Chapin III and Shaver, 1985;Shaver et al., 1986;Johnson and Caldwell, 1975) due to the unique prevailing environmental conditions.In wetland regions, water can both stimulate and limit plant productivity.Snowmelt provides a large portion of annual precipitation to Arctic regions, and the high humidity of growing season conditions limits water loss.Water tables are above or at the ground surface in moist-wet tundra ecosystems, and beneath or at the rooting level in shrub-dry tundra (Chapin III et al., 2000).Permafrost both limits percolation past the rooting depth, and provides an added input of water to plant roots throughout the growing season (Oberbauer and Dawson, 1992).In low Arctic regions, water availability is therefore not linearly associated with plant productivity (Oberbauer and Miller, 1979;Miller, 2006;Chapin III and Shaver, 1985).
In polar desert regions, surface drying can occur despite ongoing saturation of sub-surface soils, and surface drying therefore has little biological influence on plant productivity (Gold and Bliss, 1995).Water does have an indirect influence in determining Arctic vegetation species distributions due to its role in germination (Bliss, 1958); however, in PolarVPRM, Arctic tundra vegetation remains within the same allocated vegetation class (i.e., graminoid tundra, shrub tundra, or barren/wetland) throughout the model runs (< 20 years).W scale is therefore always set to 1 for regions with tundra vegetation, and is calculated according to LSWI in forested areas of North American high-latitude (NAHL).

Respiration
VPRM and PolarVPRM simulate R as a function of temperature, where R encompasses autotrophic and heterotrophic respiration.During the growing season, vegetation growth and maintenance respiration have a larger contribution to R than soil respiration does.In other words, growing season R is more heavily influenced by aboveground than belowground temperatures.Growing season R in VPRM and Po-larVPRM is therefore simulated as a function of air temper-ature.In VPRM, R is estimated year round as a piecewise linear function of air temperature, meaning that R is set to a low constant value throughout the portion of the year when air temperatures are low.
In PolarVPRM, snow season R is calculated according to soil temperature, rather than air temperature, because rates of subnivean respiration are driven primarily by soil temperature rather than air temperature (Grogan and Jonasson, 2006;Panikov et al., 2006;Sullivan et al., 2008;Morgner et al., 2010).Arctic field studies have shown that a large portion of annual carbon efflux can occur during the snow season (Aurela et al., 2004;Sullivan et al., 2008;Elberling and Brandt, 2003), and that the low thermal conductivity of an overlying snowpack (e.g., 0.06 W (m • C) −1 ; Sturm, 1992) substantially decouples Arctic soil and air temperatures, e.g., by 10-40 • C; (Zimov et al., 1993), or by 15-20 • C; (Olsson et al., 2003).Calculating snow season R according to soil temperature, rather than setting R to a low constant value like in VPRM, is therefore likely to better capture inter-annual and seasonal variability in snow season NEE, and reduce uncertainty in annual Arctic C budgets (Luus et al., 2013c).
Accuracy in estimates of Arctic R throughout the snow and growing seasons is maximized by first demarcating the snow and growing seasons according to MODIS observations of fractional snow cover area (SCA) (Appendix A), since filtered MODIS estimates of SCA agree well with in situ observations of SCA (Luus et al., 2013a).
Snow season (SCA ≥ 50 %) respiration is then calculated as a linear function of soil temperature, and growing season respiration (SCA < 50 %) is calculated as a piecewise linear function of air temperature: Calculating subnivean respiration from soil temperature and growing season respiration from air temperature decreased model errors at two high-latitude sites (Daring Lake and Ivotuk), relative to other model formulations using either soil or air temperatures alone, including the original VPRM (Luus et al., 2013a).At all three Arctic calibration sites, R 2 values were larger when respiration was regressed linearly from air temperature, than when exponential or Q 10 functions were used to describe these associations, likely because Arctic rates of respiration are low, and therefore only the low part of the exponential curve is captured.Furthermore, statistical tests with data from Ivotuk (2004Ivotuk ( -2007) ) using Akaike's information criterion (AIC) found lower AIC scores when respiration was estimated from air and soil temperatures, than when respiration was estimated from air temperature alone.These AIC scores indicate that model quality is improved by the inclusion of soil temperature, despite the concurrent increase in model complexity.

PolarVPRM parameterization by vegetation class
High-latitude vegetation is heterogeneous, resulting in large variability in NEE and its drivers by vegetation type (Humphreys and Lafleur, 2011;Elberling, 2007).Po-larVPRM therefore separates high-latitude vegetation into seven classes using a combination of the Synergistic Land Cover Product (SYNMAP) (Jung et al., 2006) and the circumpolar Arctic vegetation map (CAVM) (Walker et al., 2005) (Table 1).CAVM is available only above the northernmost tree line, whereas SYNMAP is available globally.CAVM estimates are therefore used wherever available, and SYNMAP estimates are used to classify vegetation south of the CAVM tree line.The combined CAVM-SYNMAP vegetation classification is available at 1 km × 1 km resolution, and can be upscaled when coarser resolution (≥ 1 km 2 ) Po-larVPRM outputs are desired.
Levene's test was previously applied to determine whether the CAVM-SYNMAP classes in Table 1 delineate pan-Arctic groupings with heterogeneous distributions of passive microwave-derived estimates of Arctic NEE drivers.These included cold season snow water equivalent, and growing season soil moisture, air temperature, and vegetation opacity.Findings indicated that for all passive microwave estimates examined, the CAVM-SYNMAP classes divided the pan-Arctic region into distributions with heterogeneous variances (p value < 10 −5 ), and which all displayed homoscedasticity over time (0.5 ≥ p value ≥ 0.99) (Luus et al., 2013b).The distinct distributions of the vegetation classes used, and their stability of variances over time, both indicate suitability for modeling purposes.
The NAHL spatial resolution selected for this project is 1/6 • × 1/4 • (latitude × longitude), and therefore vegetation classes are regridded accordingly.Each pixel is characterized by its fractional cover by one or more vegetation classes, and by its fractional water/glacier cover.NEE is calculated separately for each vegetation class, and total NEE for each pixel is calculated by multiplying the NEE for each vegetation class by its fractional cover.
Six separate parameters are used in the calculation of NEE across each vegetation class.All parameters are set empirically from associations found between meteorological and EC tower observations at one calibration site per vegetation class.Two parameters are used to calculate GEE (Eq.5), and four parameters are used to calculate respiration (Eq.6): -PAR 0 describes the sensitivity of photosynthetic uptake to the quantity of incoming shortwave radiation.
λ represents light-use efficiency of vegetation, and also acts as a scaling parameter.
α and β regression coefficients describe the linear association between growing season respiration and air temperature.α s and β s regression coefficients describe the linear association between soil temperature and snow season respiration.

PolarVPRM inputs
PolarVPRM remote-sensing observations of the land surface were acquired from MODIS, and meteorological observations were acquired from the North American regional reanalysis (NARR), at the native spatial and temporal resolutions summarized in Table 2.All inputs were then regridded to the PolarVPRM spatial domain (NAHL, 1/6 • ×1/4 • -latitude × longitude) using bilinear interpolation in R. MODIS was linearly interpolated to 3-hourly time steps.PolarVPRM was run at NARR's native temporal resolution (3 hourly), although it could easily be run at a finer spatiotemporal resolution with higher-resolution inputs.For a complete discussion on the changes made to model inputs, and reasons why specific meteorological products were selected, please refer to Appendix A.

Calibration and validation sites
PolarVPRM was calibrated and validated using standard meteorological observations and open-path eddy covariance measurements of NEE collected at high-latitude (HL) North American sites (Table 3; Fig. 1).All parameters except λ were set according to half-hourly EC and meteorological observations, and λ was set using observations averaged to 3hourly timescales to match the temporal resolution of Po-larVPRM.
Atqasuk (AT) and Barrow (Ba) are both located on the North Slope of Alaska, and were designated as paired calibration/validation sites representing barren/wetland vegetation in PolarVPRM.Field observations at the main Barrow EC site have indicated full flooding following snowmelt, with vegetation that consists mainly of wet sedges, moss, lichens and grasses (Oechel et al., 1995;Harazono et al., 2003).The Atqasuk EC site is located ≈ 100 km south of Barrow, and is both warmer and drier than the Barrow site.The predominant vegetation at Atqasuk is moist-wet sedge, underlain by wet, acidic soils (Kwon et al., 2006).Due to the similarity of these sites, they have previously acted as paired sites in studies of the Arctic carbon cycle (Hollister et al., 2005;Huemmrich et al., 2010).
Imnavait (Im) and Ivotuk (IV) are moist tussock tundra sites which were paired for calibration/validation, and represent the graminoid tundra vegetation class in PolarVPRM.Imnavait is located in the foothills north of the Brooks Range, in a region largely dominated by moist, acidic tussock tundra (Euskirchen et al., 2012).Ivotuk is a moist, acidic tussock tundra site located ≈ 200 km south of Atqasuk (Thompson et al., 2006;Laskowski, 2010), and Imnavait is located ≈ 200 km east of Ivotuk.The similarities in predominant vegetation, and geographical proximity, allowed these two sites to be paired.
The main Daring Lake (DL) EC site is located in a region of mixed tundra, in Canada's Northwest Territories (Lafleur and Humphreys, 2008;Humphreys and Lafleur, 2011).Observations from this site were used to calibrate the shrub tundra vegetation class.Since no other year-round EC observa-  Lafleur et al. (2012) at Canadian High Arctic sites: Lake Hazen (lh), Cape Bounty (cb), Pond Inlet (pi), Iqaluit (iq).Observations from these sites were used as model validation for the graminoid and shrub tundra classes.

Methodology
Briefly, PolarVPRM estimates of 3-hourly NEE were validated against observations from nine North American sites, and a detailed error attribution was then conducted using observations from two validation EC tower sites.Output from PolarVPRM and two existing models (CarbonTracker and FLUXNET Multi-Tree Ensemble) were then compared relative to EC observations.Changes over time (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) in PolarVPRM estimates of carbon cycling were then examined at various spatial scales across the entire high-latitude (north of 55 • N) North American terrestrial region, hereafter referred to as NAHL, for the 2001-2012 time period.

Calibration
Sites with year-round eddy covariance observations (Fig. 1) were first described using the combined CAVM and SYN-MAP classification (Table 1), and then paired.The CAVM and SYNMAP combined classification was specifically designed to allow for ecological differences resulting in varying flux drivers to be well represented, while ensuring that no category would be created that could not be parameterized using the existing eddy covariance infrastructure.If more year-round flux towers existed, then further distinctions could have been made between vegetation classes (e.g., barren/wetland).
Atqasuk, Ivotuk and Daring Lake were then selected as the calibration sites because they alone shared year-round observations during a single common year (2005).Parameter values for PolarVPRM's three Arctic tundra vegetation classes were then set using half-hourly observations from EC and meteorological towers collected at Daring Lake, (Lafleur and Humphreys, 2008;Humphreys and Lafleur, 2011), Ivotuk (Laskowski, 2010) and Atqasuk (Laskowski, 2010) (Fig. 1).In all cases, observations of NEE were filtered only to remove observations collected during instrument malfunction or when frictional velocity was low (u * < 0.2) (Goulden et al., 1996).No gap filling was carried out for any of the EC measurements, as gap filling requires the application of another model and therefore does not represent a direct measurement of CO 2 flux (Barr et al., 2004).Map of all North American calibration and validation sites and their predominant vegetation types: graminoid tundra (red), shrub tundra (cyan) or wetland (purple).Calibration sites are indicated in all caps (e.g., AT), year-round validation sites are capitalized (e.g., Im), and growing season validation sites appear in lowercase (e.g., cb).For a summary of all study site locations, please refer to Table 3. Respiration parameters were set using simple linear regressions of air temperature and nighttime NEE (α and β), or subnivean soil temperature and snow season NEE (α s and β s ).PAR 0 was set according to a non-linear least squares fit of PAR and GEE, using the nls function in R (R Development Core Team, 2011).The light-use efficiency and scaling parameter (λ) was set to be equal to the slope from a linear regression of PolarVPRM GEE vs. daytime growing season NEE, and was jointly optimized with PAR 0 .
These parameters remained unchanged for all simulations, and were applied to generate regional estimates, as well as estimates at calibration and validation sites.Parameters for vegetation classes south of the tree line were set according to the VPRM parameterizations found in Mahadevan et al. (2008).Please refer to Mahadevan et al. (2008) for a detailed description of the eddy covariance study sites used for calibration/validation purposes, the calibration approach used and results indicating good predictions of monthly NEE over forested eddy covariance calibration sites, and their respective cross-validation sites.
Following model calibration at high-latitude sites, Po-larVPRM estimates of GEE, respiration and NEE were generated for the entire North American region north of 55 • N for years 2001-2012 at a 3-hourly time step and a spatial resolution of 1/6 • × 1/4 • (latitude × longitude).The output from these simulations was used to conduct error analysis, validation, model inter-comparisons and trend analysis.

Validation
Validation consisted of examining model performance both over paired calibration/validation sites (AT, Ba, IV, Im), as well as growing season validation sites (lh, cb, pi, iq, ch).These sites capture a wide variety of vegetation types and regions of the North American Arctic, especially in light of the small total number of sites in this region with continuous year-round observations during the MODIS era (2000-).
Model evaluation consisted of examining the mean bias error (MBE) and root mean squared error (RMSE) between PolarVPRM estimates of mean 3-hourly, daily and monthly average NEE (pred i ), and EC measurements of NEE (obs i ) at matching temporal resolutions: Validation was conducted against 2005 EC observations at the three calibration sites (DL, AT, IV), and against observations from 2008 and 2001 at Im and Ba, respectively; 2008 and 2001 were selected as these were the closest years to 2005 for which year-round observations existed.After error metrics were calculated, plots comparing modeled and observed values of NEE were generated to assist in identifying biases in modeled NEE.Validation was also conducted using observations of NEE collected in July 2008 by Lafleur et al. (2012) from four Canadian Arctic sites (cb, iq, lh and pi).
Mean daily comparisons of NEE were used because of the large number of gaps found in eddy covariance observations, and the decision to not apply any gap-filling approaches to the flux data, as this would then constitute an inter-model comparison rather than a comparison of model estimates against eddy covariance observations.To ensure that these daily flux estimates would not be biased relative to model estimates in situations where more gaps in flux observations occurred at night rather than during the day, model estimates corresponding to time periods with missing observations were not included when calculating mean daily NEE.Validation therefore described the fit between model output, and in situ observations of NEE.

Error analysis
Due to the simple mathematical formulation of VPRM and PolarVPRM, uncertainties in estimates of NEE can be easily partitioned into systematic versus random errors (Lin et al., 2011).Systematic errors or biases cause model output to be offset in a specific direction, whereas random errors introduce additional and erroneous fluctuations in value.In order to better understand the deviation between PolarVPRM estimates of NEE and EC measurements of NEE, a comprehensive error analysis was completed according to the framework developed by Lin et al. (2011), which is based on a firstorder Taylor expansion.
Within this framework, errors are quantified by examining model estimates against eddy covariance observations at two year-round validation sites, Imnavait and Barrow, which had not been used in model calibration.Errors are then classified as either systematic or random.Errors are then attributed to input variables and parameters, and their total contributions to uncertainty in estimates of NEE are examined.Biases occurring due to input variables are addressed by comparing inputted shortwave radiation, soil temperature and air temperature, to NARR-derived estimates of these variables, and examining the portion of error in NEE arising from these discrepancies.Typically, model runs rely on using parameters fitted at the calibration sites (Ivotuk and Atqasuk).Biases occurring due to mis-parameterization were assessed by first fitting all parameters using EC and meteorological observations from the validation sites (Im and Ba), then comparing model NEE, generated using calibration-site parameters (IV and AT), to model NEE generated using site-specific parameters (from Im and Ba, respectively).Plotting the contribution of each component to model error then allows their relative contributions to be assessed.

Model inter-comparison
PolarVPRM estimates of NEE were compared against estimates of NEE by existing models with different formula-tions.All models were compared against EC NEE, which was upsampled from a half-hourly temporal resolution to 3hourly, daily and monthly time periods.Daily and 3-hourly averages were created using only model estimates for which concurrent EC NEE observations had been collected, in order to complete analysis without gap-filled EC NEE data.
The models selected for inter-comparison were Carbon-Tracker CT2011_oi and FLUXNET model tree ensemble (MTE).CarbonTracker (Peters et al., 2007) and FLUXNET MTE (Jung et al., 2009) were selected on the basis that both provide estimates of NEE over northern regions, and both use approaches which are different from one another and from PolarVPRM.
CarbonTracker derives estimates of CO 2 surface fluxes by analyzing atmospheric CO 2 observations using a transport model (Transport Model 5, TM5) (Krol et al., 2005) in combination with a land surface biospheric flux model (Carnegie-Ames-Stanford approach, CASA) (Potter et al., 1993) and fossil fuel inventories.One identified source of uncertainty in CarbonTracker estimates is from measurement errors or biases in CO 2 dry mole fractions (Masarie et al., 2011).FLUXNET MTE generates regional estimates of NEE by first training an ensemble of model trees using EC measurements from FLUXNET sites and inputs from the Lund-Potsdam-Jena managed Land (LPJmL) model (Bondeau et al., 2007), and then upscaling these measurements accordingly (Jung et al., 2009).
Estimates of NEE by PolarVPRM, CarbonTracker and FLUXNET MTE were spatially averaged for the entire terrestrial region north of 55 • N at a monthly resolution for each year (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).Visual examination of Fig. 3, showing monthly average NEE for each model, provided insights into differences in high-latitude carbon cycling estimated by these models.Further analysis then consisted of creating similar plots examining PolarVPRM and CarbonTracker model output at distinct time slices (Fig. 4), as both of these models are available at a 3-hourly resolution.In order to create these plots, a distinct time of day was selected for each plot, and model estimates were then spatially averaged over the North American high-latitude domain, and averaged for each month of each year.Analysis of these plots highlighted differences in the representation of diurnal patterns in highlatitude NEE.
PolarVPRM, CarbonTracker and FLUXNET MTE estimates of 3-hourly and mean monthly NEE were then compared against 3-hourly and monthly observations of NEE available at the PolarVPRM calibration and validation sites for which annual observations were available: Atqasuk, Barrow, Daring Lake, Imnavait and Ivotuk.

Trends (2001-2012)
Changes over time in the high-latitude carbon cycle were first examined by plotting total CO 2 exchange across the NAHL model domain to determine the relative contributions of res-   piration and photosynthesis (Fig. 5).Although recognized model uncertainties and the impossibility of thoroughly evaluating PolarVPRM performance across the heterogeneous model domain with current infrastructure limit confidence in estimates of the total carbon balance, examination of relative changes in CO 2 exchange over time and its drivers provide insights into responses of the high-latitude carbon cycle to recent environmental changes.
Trends over time were examined first for each year and each vegetation class, and then pixel-by-pixel across the entire model domain .The non-parametric slope estimator (Sen's slope) (Sen, 1968) was applied to each pixel to determine the trend during the years 2001-2012 of NEE, GEE and respiration.All calculations of Sen's slope values and their significance were conducted in R using the rkt package (Marchetto, 2012).
To further understand the specific influences driving these shifts, changes over time in carbon cycle variables and driver data were separately analyzed over the snow season (SS) and growing season (GS).These time periods were differentiated using MODIS MOD10A2 SCA, as a previous study indicated that remote-sensing estimates of 50 % SCA can accurately capture the timing of seasonal transitions in the low Arctic (Luus et al., 2013a).Changes over time in model inputs/outputs describing land surface characteristics are therefore described annually, as well as over the snow season (when SCA ≥ 50 %).Estimates of model variables over the growing season are limited to the portion of the year for which the land surface is snow free and for which vegetation is photosynthetically active (SCA < 50 % AND GEE < 0).This was done so that different insights could be gained into annual GEE and respiration, vs. GEE and respiration during the active growing season.The Sen's slope estimates of median changes in carbon cycle and land surface variables over time (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) were then reported for each pixel in the model domain corresponding to significant (p value < 0.05) change.

Validation
PolarVPRM 3-hourly estimates of NEE showed excellent agreement with 3-hourly averaged NEE measured at four Canadian Arctic sites during the 2008 growing season (Lafleur et al., 2012).The mean RMSE across all sites was 0.79 µmol CO 2 m −2 s −1 , and was largest at Iqaluit (1.01 µmol CO 2 m −2 s −1 ) and smallest at Cape Bounty (0.66 µmol CO 2 m −2 s −1 ) (Table 4).Although these errors are not negligible, the values are quite low considering the large distances between calibration sites and growing season validation sites (Fig. 1), as well as the large inter-site environmental differences (Lafleur et al., 2012).
Biases in NEE at these sites arose mainly due to biases in NARR shortwave radiation and air temperature.At Lake Hazen, NARR air temperatures were 10 • C colder on average than observed air temperatures, leading both photosynthesis and respiration to be underestimated.Furthermore, EVI remained low at Lake Hazen (< 0.1) and Iqaluit (< 0.2) throughout July 2008.Photosynthesis at these sites, and total net C uptake in July 2008, was therefore underes-q q q q q q q q q q q q q q q q q q q q q q q q Respiration −1 x GEE timated at these sites (lhMBE = 0.39, iq MBE = 0.15).At Cape Bounty and Pond Inlet, NARR shortwave radiation was substantially larger than observed, such that cb measured PAR ≈ 1.3 × cb NARR SW radiation, and iq measured PAR ≈ 1.6 × iq NARR SW radiation.These biases in shortwave radiation cause overall rates of photosynthesis and net C uptake at these sites to be overestimated (cb MBE = −0.33,pi MBE = −0.44).Despite the small biases introduced from meteorological inputs, PolarVPRM growing season observations agree relatively well with observations from remote Canadian Arctic sites (mean RMSE = 0.79).

Error analysis
Comparisons of PolarVPRM NEE against observations of NEE at year-round calibration and validation sites indicated that PolarVPRM tended to slightly underestimate net C uptake by vegetation at validation sites (Barrow and Imnavait).
The underlying reasons for the bias in net C uptake at yearround validation were addressed through a comprehensive error analysis.
The cumulative monthly bias in PolarVPRM, expressed as the difference between PolarVPRM modeled NEE and observed NEE (model−observed), is indicated for two validation sites, Barrow and Imnavait in Fig. 2. In this figure, a positive bias in NEE indicates that either respiration was overestimated, or that photosynthetic uptake by vegetation was underestimated.Therefore, according to the sign con-vention used for NEE, a positive bias in respiration indicates that PolarVPRM overestimated respiration, whereas a positive bias in GEE indicates that PolarVPRM underestimated photosynthetic uptake by vegetation.
The main sources of error in PolarVPRM arose from biases in how the associations between PAR and GEE, and between λ and GEE, were parameterized (Fig. 2).NARR air temperature and soil temperature agree closely with observed soil and air temperatures, meaning there is very little bias in estimates of respiration.Small errors in NARR shortwave radiation relative to observations caused small (< 1 g C m −2 ) cumulative biases in net C exchange at these validation sites.
At Ba and Im, the amount of carbon taken up through photosynthesis is underestimated due to biases in λ, the lightuse efficiency and scaling parameter (Eq.5).At the calibration sites, AT and IV, λ values of 0.15 and 0.04 were identified as being optimal values for barren/wetland regions and graminoid tundra sites, respectively.When optimal λ were instead calculated using EC NEE from validation sites (Ba A recent study by Dietze et al. (2014) likewise identified mis-parameterization of light-use efficiency at low-light levels to play a central role in biospheric model uncertainties across HL regions of North America.The study indicated that a likely source for mis-parameterization is due to greater variance in this parameter across high-latitude sites, even when these contain similar biota.Although PolarVPRM calibration and validation sites were paired on the basis of having similar physical and biological characteristics, important differences appear to exist in the drivers of NEE at these sites, which are not captured by PolarVPRM.In order to identify whether the calibration or validation LUE values are more representative across HL tundra sites, a larger network of HL EC towers and measurements of light-response curves using leaf-level observations of gas exchange (Bernacchi et al., 2013) would be required.
A portion of the error observed at Barrow and Imnavait arises due to biases which are not considered in the error analysis framework.The gap between the sum of all biases and the total error is evident at both Barrow and Im-navait (Fig. 2).The error analysis approach used is based on a first-order Taylor expansion, and therefore does not consider second-order effects (Lin et al., 2011).As PolarVPRM has a simple model structure, it is likely that a portion of the bias arises due to incomplete characterization of the processes which drive carbon cycling.Future work will consist of further improving the accuracy of PolarVPRM estimates.Spatial heterogeneity in rates of respiration and CO 2 release from permafrost could be better represented by including satellite-derived maps of permafrost information (e.g., Heim et al., 2011).Estimates of PolarVPRM GEE could also be improved in the future once capabilities exist to accurately estimate vegetation fluorescence or photosynthetic stress from the photochemical reflectance index (Grace et al., 2007;Hilker et al., 2008) across vast Arctic regions.Overall, although mis-parameterization of λ at the validation sites accounts for the most significant portion of error at Barrow and Imnavait, a hidden bias in NEE exists that compensates for the bias in λ, resulting in smaller net biases in NEE than would be expected through error decomposition.

Model inter-comparison
PolarVPRM shows closer agreement with EC NEE from five Arctic sites, than FLUXNET MTE shows against the same five sites (Table 5), indicating that PolarVPRM provides an improved data-driven approach for estimating regional-scale Arctic NEE.When 3-hourly, daily and monthly averages of All statistically significant (p value < 0.05) changes over time in carbon cycle variables and driver data are shown for the growing season (GS, when SCA < 50 % AND GEE < 0).Please note that as the growing season includes only the period of time for which vegetation is productive at any pixel (GEE < 0), periods of time for which air temperature extremes or drought hinder photosynthesis are not included.
The influences of rising EVI and air temperatures on increasing Arctic rates of photosynthesis are therefore made clear, whereas plots of annual GEE (Fig. 7) emphasize reductions in photosynthetic uptake of C by forest vegetation.
from five sites at the same timescales, PolarVPRM had the lowest mean RMSEs for all timescales, and lower MBEs at monthly timescales, but larger MBEs at daily and 3-hourly timescales.PolarVPRM therefore provides estimates of NEE which show similar or improved realism relative to EC NEE, using a simpler framework than CarbonTracker.
Comparisons of mean monthly NEE by PolarVPRM, Carbon-Tracker and FLUXNET model tree ensemble over NAHL indicated that both CarbonTracker and PolarVPRM estimated very low rates of mid-winter respiration, whereas the FLUXNET MTE showed greater rates of mid-winter respiration (Fig. 3).FLUXNET MTE also estimates greater photosynthetic uptake of carbon by vegetation than the other two models.The seasonal cycle and inter-annual variability displayed by PolarVPRM and CarbonTracker are very similar (Fig. 3).
Relative to CarbonTracker, PolarVPRM estimates less carbon to be taken up by vegetation photosynthetically at midday (18:00 UST), and estimates less respiration to occur in the middle of the night during the growing season (06:00 UST) (Fig. 4).However, when 3-hourly estimates of NEE by PolarVPRM and CarbonTracker are compared against 3-hourly EC measurements of NEE, PolarVPRM is found to have lower RMSE values than CarbonTracker (Table 5).CarbonTracker has a lower summed MBE than Po-larVPRM because PolarVPRM more substantially underestimates peak growing season GEE at Daring Lake than Car-bonTracker does, and this bias is disproportionately large as EC NEE is only available at Daring Lake during the growing season (unlike other sites).If 3-hourly and daily MBEs were instead summed over the remaining sites, PolarVPRM would have lower MBEs than CarbonTracker.
PolarVPRM may have lower RMSEs than CarbonTracker because CarbonTracker estimates less photosynthesis to occur at 00:00 and 12:00 UST than estimated by PolarVPRM.Arctic regions receive sunlight continuously through the mid-summer season, and Arctic tundra vegetation has been observed to conduct photosynthesis continuously throughout mid-summer, despite low-light levels and low temperatures found at solar midnight (Tieszen, 1973;Patankar et al., 2013).PolarVPRM's ability to simulate photosynthesis according to actual light availability at high latitudes therefore allows estimates of NEE to agree closely with averaged EC NEE, especially when considered at monthly intervals.In summary, PolarVPRM provides a simple approach for generating reliable estimates of NEE.

Regional trends (2001-2012)
Regional trends are examined here for the purpose of examining general tendencies in inter-annual variability, as uncertainties in regional-scale estimates of North American NEE limit the capacity to accurately quantify high-latitude regional-scale C balances.Model estimates of the mean annual carbon balance of the North American region north of 55 • N indicate that this region may have been a carbon source from 2001 to 2004 and 2010 to 2012, and a weak carbon sink between 2005 and 2009 (Fig. 5).Between 2001 and 2009, respiration diminished, and then rose again from 2010 onward.After 2004, GEE became less negative over time, indicating less CO 2 uptake per year.Respiration changes more than GEE over time, and appears to have a larger role in determining changes over time in NAHL net C uptake than photosynthesis.
Several general tendencies appear when examining monthly average NEE over time (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) from each PolarVPRM vegetation class (Fig. 6).In 2001-2005, Po-larVPRM estimates that high-latitude North America was a carbon source (Fig. 5) as the source strength of tundra regions exceeded the sink strength of forested regions.Over time (2005)(2006)(2007)(2008)(2009)(2010), forested regions took up less carbon through photosynthesis year after year (Fig. 6).Graminoid tundra and barren regions continue to function as net carbon sources, and forested regions continue to function as net carbon sinks.Shrub tundra regions appear to be shifting towards becoming carbon sinks because of an increase in the amount of C taken up through photosynthesis.The decrease in respiration observed in Fig. 5 is due largely to a small decrease in snow season respiration.Although this decrease is minor on a monthly level, its cumulative impact on the carbon balance is substantial.Forests have stronger fluxes than tundra vegetation, and therefore have a greater relative contribution to the North American high-latitude carbon cycle than tundra regions, leading to a net increase in CO 2 efflux from 2007 onward (Fig. 5).

Per-pixel trends (2001-2012)
Previous studies have described the influence of warming air temperatures on inducing increased rates of net carbon uptake by vegetation near the shrub and tree lines (Hinzman et al., 2005;Tape et al., 2006), and on increasing rates of CO 2 efflux (Schuur et al., 2009;Tarnocai, 2006).Remotesensing studies have found trends towards increased growing season length (Zeng et al., 2011), increased normalized difference vegetation index (NDVI) over tundra regions due to warming (Stow et al., 2004), and diminished NDVI over boreal regions due to reduced rates of photosynthesis (Verbyla, 2008).Since PolarVPRM is driven by remote-sensing observations, the effects of these environmental changes for HL carbon cycling can be examined by analyzing trends in Po-larVPRM output.
Trends over time in carbon cycle variables were examined for each pixel in NAHL individually using the nonparametric Theil-Sen estimator (Sen's slope).Initial analyses were conducted according to mean annual values of net carbon uptake, NEE, respiration and GEE (Fig. 7).Visual examination of these plots indicated a net increase in carbon efflux from high-latitude regions, focused mainly in forested regions (Fig. 7a).The observed increase in annual carbon efflux from forested regions over time arises mostly due to a decrease over time in the photosynthetic uptake of carbon (represented by an increase in GEE) (Fig. 7c).Although photosynthetic uptake in parts of northern Alaska and the Yukon increased over time, greater uptake was outweighed by the declines observed over forested regions.The net change in GEE is therefore mainly indicative of diminished sink strength over time.Effluxes of CO 2 from Arctic tundra regions increased over time due to greater rates of respiration (Fig. 7b).Overall, this results in a slight trend toward less net CO 2 uptake across the entire model domain (Fig. 7a), especially from 2005 onward.
The amount of carbon taken up by vegetation through photosynthesis increased over tundra regions, and declined sharply over forest regions (Fig. 7c).When considering trends only over the active growing season (when GEE < 0), there was a slight increase in the amount of carbon taken up by North American vegetation during the growing season (Fig. 8a).This discrepancy is due to the inability of model vegetation to conduct photosynthesis when temperatures rise above the maximum air temperatures permitted for photo- synthesis, and due to increased drought stress in warm conditions.As air temperatures warm above the physiologically optimal temperatures, and drought stress increases, the capacity for photosynthesis diminishes strongly.
Positive air temperature anomalies and increased drought stress during the growing season therefore limit the total amount of carbon taken up by forest vegetation, while generally, rising EVI (Fig. 8c) and air temperatures (Fig. 8d) increase photosynthetic activity whenever temperatures are not excessively hot.It is also interesting to note that a further consequence of rising air temperatures is a concurrent rise in growing season rates of respiration (Fig. 8b), which seems to partially counteract the increase in photosynthesis observed across the model domain (Fig. 8a).
Cooling trends in snow season NARR T air resulted in colder T soil over time south of the tree line, resulting in diminished respiration, but less change over time was observed in soil and air temperatures north of the tree line (Fig. 9c  and d).Recent studies have indicated cooling trends in winter air temperatures over high-latitude regions of North America could be due to trends in the Madden-Julian oscillation (Yoo et al., 2011), or due to deepened Eurasian snow depth (Cohen et al., 2012).As snow season length also diminishes over time (Fig. 9b), it could be expected that subnivean effluxes of CO 2 would contribute less carbon annually to the atmosphere over time.Only a small decline over time was observed in snow season respiration over forested regions (Fig. 9a).Conversely, diminished snow season length could contribute to the observed rises in growing season respiration (Fig. 8b).The initial decline, and later rise, in respiration are therefore likely to occur due to counteracting trends over time in snow and growing season respiration.
The boreal forest appears to have a dominant role in determining fluxes of CO 2 over North American latitudes north of 55 • N.Although photosynthetic uptake in tundra regions increases over time, this is largely outweighed by concurrent rises in respiration due to warming air temperatures.Forest regions are also capable of greater rates of photosynthesis during the active growing season (GEE < 0), but carbon uptake is limited due to drought and temperature stress.As a result, reductions over time occur in the amount of carbon taken up by vegetation.Furthermore, although subnivean effluxes of CO 2 diminish over time due to shortened snow seasons and diminished snow season soil temperatures, annual rates of respiration increase over time.PolarVPRM simulations indicate that it is likely that North American HL regions have recently been emitting more CO 2 into the atmosphere in response to warming air temperatures.

Conclusions
PolarVPRM provides a remote-sensing-based approach for generating high-resolution estimates of NEE using a parsimonious model approach that is specifically adapted for high-latitude regions.PolarVPRM adequately simulates high-latitude NEE by using Arctic-specific vegetation classes, and calculating snow and growing season R separately using soil and air temperatures, respectively.When PolarVPRM, CarbonTracker and FLUXNET MTE were all examined against EC NEE averaged over 3-hourly, daily and monthly timescales, the smallest RMSE values for each timescale were found in comparisons of PolarVPRM NEE to EC NEE, indicating reasonable model performance.
Initially, high-latitude regions increased their net uptake of carbon over time (2001)(2002)(2003)(2004)(2005) due to an increase in rates of photosynthesis by Arctic vegetation.Subsequently, net carbon efflux from high-latitude regions increased (2011)(2012) due to declines in photosynthesis over boreal regions in response to temperature and drought stress.Overall, PolarVPRM indicates that warmer air temperatures are enabling Arctic vegetation to take up more carbon photosynthetically, while simultaneously increasing high-latitude rates of respiration, and diminishing photosynthetic uptake of carbon by boreal vegetation.

Figure 1 .
Figure1.Map of all North American calibration and validation sites and their predominant vegetation types: graminoid tundra (red), shrub tundra (cyan) or wetland (purple).Calibration sites are indicated in all caps (e.g., AT), year-round validation sites are capitalized (e.g., Im), and growing season validation sites appear in lowercase (e.g., cb).For a summary of all study site locations, please refer to Table3.

Figure 2 .
Figure 2. Cumulative monthly bias in PolarVPRM estimates of net C exchange at Barrow (left) and Imnavait (right), relative to eddy covariance observations at these sites.Errors in GEE are designated as being due to the associations between GEE and downward shortwave radiation (SW_GEE), GEE and light-use efficiency (LUE_GEE), and of the parameter describing the association between GEE and PAR0 (PAR0_GEE).Shaded areas surrounding (PAR0_GEE) and (LUE_GEE) represent the range of biases possible from the determination of PAR0 and λ from eddy covariance observations.Total biases in temperature and GEE (T_GEE), and between temperature and respiration (T_R) are also described, along with the total biases in respiration (R_all) and NEE.Comparisons are shown for the range of months for which eddy covariance observations were acquired at Barrow in 2001 (January-September), and at Imnavait in 2008 (February-October).

Figure 4 .
Figure 4. Spatially averaged North American high-latitude NEE from PolarVPRM (blue) and CarbonTracker (red).All estimates are averaged monthly at four distinct times of day, shown in universal time (UST) rather than according to local time zones.

Figure 5 .
Figure 5. Relative contributions of PolarVPRM respiration and photosynthesis (GEE, plotted here as −1×GEE) to inter-annual variability in the net C balance of the North American terrestrial region north of 55 • N (NAHL).

Figure 6 .
Figure 6.Average monthly NEE (in µmol CO 2 m −2 s −1 ) for the entire North American region north of 55 • N (Mean), as well as for seven vegetation classes within this region: deciduous forest (DECDS), evergreen forest (EVGRN), mixed forest (MIXED), shrub (SHRUB), graminoid tundra (GRMTD), shrub tundra (SRBTD) and barren/wetland (BARRN).The percentage of the model domain covered by each vegetation class is indicated in the title of each subplot.

Figure 7 .
Figure 7. Sen's slope values, indicating the median change (2001-2012) in PolarVPRM estimates of mean annual NEE (a), respiration (b)and GEE (c).In (a), changes over time in NEE are indicated in bothµmol CO 2 m −2 s −1 and tC ha −1 using the same color scheme.All Sen's slope values shown correspond to p values < 0.05.Pixels with > 50 % fractional water content are indicated in grey.Please note that the negative sign convention in GEE has been maintained, meaning that a positive trend in GEE corresponds to diminished uptake of carbon through photosynthesis.

Figure 8 .
Figure 8. Sen's slope of median change (2001-2012) in PolarVPRM estimates of growing season carbon cycle variables, and driver data.All statistically significant (p value < 0.05) changes over time in carbon cycle variables and driver data are shown for the growing season (GS, when SCA < 50 % AND GEE < 0).Please note that as the growing season includes only the period of time for which vegetation is productive at any pixel (GEE < 0), periods of time for which air temperature extremes or drought hinder photosynthesis are not included.The influences of rising EVI and air temperatures on increasing Arctic rates of photosynthesis are therefore made clear, whereas plots of annual GEE (Fig.7) emphasize reductions in photosynthetic uptake of C by forest vegetation.

Figure 9 .
Figure 9. Sen's slope of median change (2001-2012) in PolarVPRM estimates of snow season (SS, when SCA ≥ 50 %) carbon cycle variables, and driver data.Values are only shown for locations with a statistically significant change over time (p value < 0.05).

Table 1 .
Jung et al. (2006)on classes, created by combining and aggregating CAVM and SYNMAP vegetation classes.SYNMAP tree classes are described according to leaf type (broad, needle or mixed) followed by leaf longevity (evergreen, deciduous, or mixed), as inJung et al. (2006).

Table 2 .
Summary of meteorological and land surface remote-sensing inputs to PolarVPRM.

Table 3 .
Calibration and validation sites for each vegetation type, years of data used, and locations.Calibration sites are in bold, and validation sites at which cold season EC observations were unavailable are italicized.
tions were available at NAHL sites designated as shrub tundra from 2001 to 2012, validation for this class consisted of characterizing model performance over years which were not used for validation, and by describing the performance of this parameterization in describing NEE at Ivotuk.Meteorological and EC observations were collected during a portion of the 2008 growing season by

Table 4 .
MBE and RMSE (in µmol CO 2 m −2 s −1 ) from comparisons of 3-hourly PolarVPRM NEE, to July 2008 eddy covariance observations of NEE at four Canadian Arctic validation sites.

Table 5 .
Error statistics (RMSE and MBE; in µmol CO 2 m −2 s −1 ) found through the comparison of monthly, daily and 3-hourly averaged estimates of NEE from PolarVPRM, CarbonTracker and FLUXNET model tree ensemble relative to observations of 3-hourly, daily and monthly averages of NEE from Atqasuk (AT), Barrow (Ba), Daring Lake (DL), Imnavait (Im) and Ivotuk (IV).