Modeling the influence of snow cover on low Arctic net ecosystem exchange

The Arctic net ecosystem exchange (NEE) of CO2 between the land surface and the atmosphere is influenced by the timing of snow onset and melt. The objective of this study was to examine whether uncertainty in model estimates of NEE could be reduced by representing the influence of snow on NEE using remote sensing observations of snow cover area (SCA). Observations of NEE and time-lapse images of SCA were collected over four locations at a low Arctic site (Daring Lake, NWT) in May–June 2010. Analysis of these observations indicated that SCA influences NEE, and that good agreement exists between SCA derived from time-lapse images, Landsat and MODIS. MODIS SCA was therefore incorporated into the vegetation photosynthesis respiration model (VPRM). VPRM was calibrated using observations collected in 2005 at Daring Lake. Estimates of NEE were then generated over Daring Lake and Ivotuk, Alaska (2004–2007) using VPRM formulations with and without explicit representations of the influence of SCA on respiration and/or photosynthesis. Model performance was assessed by comparing VPRM output against unfilled eddy covariance observations from Daring Lake and Ivotuk (2004–2007). The uncertainty in VPRM estimates of NEE was reduced when respiration was estimated as a function of air temperature when SCA ≤ 50% and as a function of soil temperature when SCA > 50%.


Introduction
In low Arctic regions, the initial onset and final melt of snow mark important transitions in net ecosystem exchange (NEE) (Olsson et al 2003, Grogan et al 2004, Bokhorst et al 2010, Buckeridge and Grogan 2010, where NEE is defined Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. as the net biospheric flux of CO 2 into and out of the land surface (Lovett et al 2006). NEE can be described as the sum of photosynthetic uptake by vegetation (GEE, or gross ecosystem exchange) and ecosystem respiration (R): NEE = −GEE + R. According to the sign convention used in this study, uptake of CO 2 from the atmosphere is represented as negative NEE, and release of CO 2 into the atmosphere is shown as positive NEE.
Photosynthetic uptake by vegetation is maximized during the growing season, when above-freezing air temperatures (T air ) and sunny conditions support plant growth. During snow onset in autumn, the land surface is cooling (Zhang et al 1996), light availability is limited, and most plants are in senescence, resulting in diminished rates of photosynthesis and respiration (Billings and Mooney 1968, Carstairs and Oechel 1978,Öquist and Huner 2003, Olsson et al 2003, Euskirchen et al 2012. Snowpack development decouples T air and soil temperature (T soil ) (Bonan 2002), allowing subnivean respiration to persist at low T air (Zimov et al 1993, Olsson et al 2003. Snow melt in the Arctic normally takes place within a month of the solstice. Light availability is therefore high, and melt is accompanied by warmer T air , soil thaw and greater availability of nutrients. As a result, snowmelt is accompanied by rapid increases in rates of respiration (Zimov et al 1996, Mikan et al 2002, Oelbermann et al 2008, Elberling et al 2008. Although the length of time between snowmelt and green-up varies by species (Humphreys and Lafleur 2011), the timing of snowmelt at a site influences the timing of photosynthetic uptake by vegetation (Morgner et al 2010, Buckeridge and Grogan 2010).
Landscape rates of photosynthesis and respiration during snow onset and snow melt are also influenced by the proportion of the land surface which is snow covered at any given point in time (snow cover area, or SCA). Comparisons of NEE at plots with differing quantities of snow have found diminished rates of photosynthesis and respiration during snow onset/melt at plots with greater SCA (e.g. Buckeridge andGrogan 2010 andMorgner et al 2010). Representing SCA in biospheric carbon flux models could therefore allow the snow, transitional and growing seasons to be clearly delimited. Model estimates of NEE during these time periods could then be generated by simulating the differing seasonal drivers of NEE for each period. Hence, uncertainty in model estimates of Arctic NEE might be reduced by explicitly representing the influence of SCA on NEE.
Although  (Dozier 1989, Rosenthal andDozier 1996) and MODIS (moderate resolution imaging spectroradiometer) (Hall et al 2002, Hall and Riggs 2007, Riggs and Hall 2011. Remote sensing observations of SCA could be applied to represent the influence of snow on decoupling soil and air temperatures (Olsson et al 2003, Bonan 2002, and causing snow season respiration to be driven by soil temperature rather than air temperature. Similarly, rates of photosynthesis are greatly limited by the presence of an overlying snowpack (Tieszen 1974), since snow scatters incoming solar radiation and therefore limits the amount of light that can reach subnivean vegetation (Warren 1982, Zhou et al 2003. Estimates of SCA from satellites such as Landsat and MODIS could therefore be incorporated into models of biospheric CO 2 fluxes, and applied to improve the representation of snow season (SCA > 0%) respiration and/or photosynthesis.
Modeling the influence of snow on NEE using a remote sensing approach offers several advantages over a process-based approach. Incorporating remote sensing observations of snow characteristics in a model of NEE limits the propagation of meteorological biases into estimates of snow characteristics, and allows spatial variability in snow distributions to be captured without having to simulate the many land surface influences on snow accumulation and melt (i.e. topography, aerodynamic roughness). In order for SCA to be included in models of NEE, it is important to first address the challenges which have previously prevented their inclusion in models of NEE. Specifically, uncertainty has existed about the accuracy of remote sensing estimates of SCA at high-latitude sites, whether the influence of snow on NEE could be detected using SCA, and whether uncertainty in model estimates of NEE could be reduced by incorporating SCA. The research provided here therefore systematically addresses these three challenges.
The central goal of this research was to explore the potential for incorporating SCA and its effects on biospheric carbon fluxes into the vegetation photosynthesis respiration model (VPRM), a diagnostic, remote sensing assimilation scheme designed to provide regional estimates of NEE (Mahadevan et al 2008). The specific objectives were to examine the feasibility of assimilating MODIS SCA into a model of biospheric carbon fluxes, and to examine whether uncertainty in VPRM estimates of NEE could subsequently be reduced by simulating the influences of SCA on photosynthesis and/or respiration.

Study sites
The Daring Lake (DL) site is located in the southern portion of the Northwest Territories at 64 • 52N, 111 • 34 W, ≈200 km northeast of Yellowknife. Daring Lake receives an average of 200-300 mm in precipitation annually, and has a mean annual temperature of −12.5 to −10.5 • C (Lafleur and Humphreys 2008). At this study site, four time-lapse cameras were automated to capture thrice daily images of the land surface in May-June 2010 from towers overlooking mixed tundra (MT), fen (FN), low shrub mixed tundra (LK), and tall shrub (SB) vegetation. Observations of NEE were simultaneously acquired from eddy covariance towers at MT and FN in May-June 2010. Measurements of NEE and meteorological variables have been collected at the Daring Lake MT site since 2004. The MT site is underlain by sand to loamy sand textured soil, and is composed of shrub tussock tundra and mesic heath. The FN site is a wet sedge meadow with 40 to 70 cm of peat soil overlying silt loam textured soil. The dominant vegetation is sedges, with minor amounts of dwarf birch and a moss understory (Humphreys andLafleur 2011, Lafleur andHumphreys 2008).
In order to select a paired validation site for Daring Lake mixed tundra (MT), all North American sites with eddy covariance and meteorological observations (2004)(2005)(2006)(2007) were considered. Of these, Ivotuk (IV) was selected as the validation site because it is the most similar to Daring Lake MT in terms of vegetation, precipitation and temperature. Ivotuk is an AmeriFlux site located on the north Slope of Alaska at 68 • 29N, 155 • 44 , ≈300 km south of Barrow. The average temperature and liquid precipitation at Ivotuk have been reported within the −8.9 to −14.6 • C and 123-221 mm ranges by Laskowski (2010). The Ivotuk site has been classified as a moist, acidic, tussock tundra site dominated by Eriophorum vaginatum and containing shrubs, mosses and lichen (Thompson et al 2006, Laskowski 2010. All eddy covariance observations of NEE were filtered to remove periods of time with low frictional velocity and time periods with sensor malfunction (e.g. when windows on the open-path infrared gas analyzer were obscured). No gap-filling was performed for observations from either site, and no additional data points were removed. The same non-gap-filled datasets were used throughout the entire study.

Calculating snow cover area
Snow cover area was estimated from visible and/or infrared observations available at three resolutions, from three different sources: time-lapse camera (<10 m), Landsat (30 m) and MODIS (500 m). Although images from both 1 May-30 June 2010 and 30 August-7 December 2010 were examined, we focus the analysis of snow-NEE associations on the time period when camera images of SCA at DL MT were acquired simultaneously with meteorological and eddy covariance observations: 1 May-30 June 2010. All camera images acquired during this two month time period were individually classified in ENVI/IDL using a combination of supervised parallelepiped and unsupervised isodata classifiers in ENVI. Unsupervised isodata classifiers were first applied. All images were then visually assessed to determine how well the isodata classifications captured SCA. For images that were not well classified using the isodata approach, a supervised parallelepiped classification was instead applied. Figure 1 shows a selection of time-lapse images from the four locations at Daring Lake at the start, middle and end of snow melt in 2010. From these analyses, a percentage of fractional snow cover was calculated three times per day.
Landsat images were collected over the 2004-2007 and 2010 periods in May-June. These images were classified in terms of snow presence/absence using the normalized difference snow index (NDSI). The NDSI has been used to identify snow on the basis that snow reflects less middle-infrared (r MIDIR , 1.55-1.75 µm) than visible (r GREEN , 0.52-0.60 µm) radiation: Regions with NDSI > 0.4 and r 0.76−0.90 µm > 11% were classified as snow covered, whereas other regions were classified as non-snow covered, as according to Hall et al (1995). The MODIS 10A1 fractional SCA product was used as the 500 m estimate of SCA. As MOD10A1 and Landsat are available on 7-9 day intervals, basic linear interpolation was used to generate three-hourly estimates of snow cover.
In VPRM, NEE has been calculated according to the sum of GEE and respiration: NEE = −GEE + RESP, and respiration has been calculated using a piecewise approach. When air T air has been warmer than a threshold temperature (T low ), respiration has been calculated as a linear function of air temperature (RESP = αT air + β). When T air < T low , respiration has been set to a low, baseline value independent of temperature. T low , α, and β have been calculated according to the relationships found between NEE and T air using tower observations. In this original formulation, cold season respiration has been assumed to remain at a constant rate throughout the time period when T air < T low , regardless of fluctuations over time in snowpack properties or subnivean temperatures.
VPRM driver data was composed entirely of remotesensing-based estimates of land surface characteristics and meteorology by MODIS and North American regional reanalysis (NARR). PAR, T air , T soil , and T scale were acquired from NARR downward shortwave radiation, 0-10 cm T soil and 2 m T air datasets. NARR estimates were generated by a model which assimilated a variety of meteorological observations (Mesinger et al 2006). At Daring Lake (2004)(2005)(2006)(2007)2010) and Ivotuk (2004Ivotuk ( -2007, NARR estimates had good agreement (R 2 > 0.8) with all meteorological tower observations of air temperature, soil temperature at 5 cm, and PAR (shortwave radiation). Preliminary investigations indicated that NARR estimates of these variables had better agreement with field observations than those generated by other leading approaches. W scale , P scale and FAPAR PAV were derived from MODIS. MODIS offers established, moderate resolution, estimates of land surface characteristics from visible and infrared observations. MODIS observations have been used extensively, including in previous versions of VPRM (Mahadevan et al 2008). Both NARR and MODIS inputs therefore represent the best available driver datasets for VPRM at high latitudes.
In this study, daily estimates of NEE in 2004-2007 at Daring Lake and Ivotuk were generated using six model formulations based on VPRM, some of which represent the influence of snow on photosynthesis (GEE s ) and/or respiration (RESP s ) according to MODIS SCA (0-100%), and some of which do not. The RESP s & GEE 0 model formulation, which calculates snow season respiration according to NARR T soil and growing season respiration according to NARR T air , was found to have the best agreement with eddy covariance observations. To determine whether this benefit arose from the calculation of respiration according to T soil rather than T air alone (and not due to the use of SCA), two model formulations that calculate respiration year-round as a linear (RESP Tsoil linear) or piecewise linear (RESP Tsoil pwl) function of soil temperature were also evaluated. The α, β, T low , α s , and β s , T Tsoil-low , α Tsoil , and β Tsoil parameters were all individually calibrated to eddy covariance observations collected at Daring Lake in 2005. In all of the model formulations below, GEE 0 refers to GEE as calculated in equation (2), and NEE = −GEE + RESP.
VPRM performance with these six model formulations was assessed both qualitatively, and statistically. The error , or difference between predicted (VPRM) and observed (non-gap-filled eddy covariance) daily average values of NEE ( i = pred i −obs i ), was evaluated using two metrics: the mean absolute error (MAE) and root mean squared error (RMSE) (Willmott and Matsuura 2005):

Landsat and MODIS estimates of local snow cover
Time-lapse camera observations of fractional SCA agreed well with linearly interpolated Landsat NDSI and MODIS observations of whether the pixel containing each camera was snow covered or snow-free (figure 2). Time-lapse camera observations and Landsat derived estimates of SCA showed snow depletion to occur within a seven day time period, with no substantial melts before snow depletion and no snow accumulation following depletion. To assess the agreement between time-lapse, Landsat, and MODIS derived estimates of the timing of snowmelt, MODIS and Landsat SCA were linearly interpolated between ≈weekly acquisitions. Comparisons indicated a slight advance or delay of one to five days in Landsat and MODIS estimates of depletion relative to time-lapse camera observations 1. A slight (<7 day) discrepancy is acceptable considering the ≈weekly temporal resolution of Landsat and MOD10A1 SCA. Overall, these results indicate good agreement between ground-based and remote sensing observations of SCA.

Associations between NEE and SCA
Preliminary investigation of SCA, T soil , T air and NEE over time at the Daring Lake MT site indicated that SCA may have several important effects on NEE (figure 3). Although the low thermal conductivity of snow enables T soil to remain warmer than T air throughout midwinter (figure 4), in the time period immediately preceding snow melt, air temperatures are warmer than soil temperatures. As T air rises, T soil slowly follows to reach a temperature of 0 • C, at which point snow melt begins. A decline in SCA during snowmelt initially increases the rate of respiration, causing an increase in CO 2 efflux. Once snow melt is complete, vegetation begins to green-up, leading to an increase in the rate of photosynthetic uptake of CO 2 . Air and soil temperatures also become more closely synchronized following snow melt. In summary, these findings suggest that the timing of snowmelt coincides closely with an increase in respiration, and subsequent increase in photosynthesis, at Daring Lake MT. It therefore appears feasible that remote sensing observations of SCA could be incorporated into VPRM to help describe the influences of SCA on respiration and/or photosynthesis.

VPRM estimates of NEE with and without MOD10A1 SCA
All VPRM formulations which estimated growing season respiration according to T air generated reasonable estimates of NEE over Ivotuk (IV) and Daring Lake MT (DL) in years 2004-2007 (MAE = 0.2-0.5 µmol m −2 s −1 and RMSE = 0.6-1.8 µmol m −2 s −1 ). MAEs were similar between DL and IV. RMSE values tend to be greater at IV than DL because a larger portion of observations at IV were collected in midwinter, a time of year in which uncertainty in observations of NEE is greatest (Amiro 2010).
Both the GEE s and GEE 0 formulations accurately represented the lack of photosynthesis during the long Arctic  midwinter. Snow melt was accompanied by a slight increase in photosynthesis (figure 3), which was best simulated by GEE 0 models. An increase in photosynthetic uptake at the end of the snow season was implicitly simulated by GEE 0 according to rising air temperatures, increased sunlight, and an eventual slight increase in EVI. Simulating further reductions in photosynthesis when 0% < SCA < 100% caused VPRM to underestimate GEE. The RMSEs and MAEs from RESP 0 & GEE s and RESP s & GEE s therefore exceeded those from the formulations which did not further suppress GEE as a linear function of SCA (RESP 0 & GEE 0 and RESP s & GEE 0 ).
It is possible that GEE 0 outperformed GEE s because, unlike GEE s , GEE 0 does not explicitly specify that no subnivean photosynthesis can take place. Indeed, subnivean photosynthesis has been observed to occur prior to snowmelt over a sub-Arctic moss-dominated heath (Larsen et al 2007). Similarly, when T air ≈ 0 • and snowcover is thin or patchy, the following Arctic species have been observed to conduct photosynthesis at <25% of the rate observed in the peak of the growing season: Carex aquatilis, Dupontia fisheri, Eriophorum angustifolium (Tieszen 1974); Eriophorum vaginatum, Ledum palustre, Vaccinium vitisidaea and Cassiope tetragona (Starr and Oberbauer 2003). A central reason why subnivean photosynthesis is a very gradual process is because snowpacks act as an optical scattering medium (Warren 1982); <5% of incoming solar radiation has been observed to penetrate >10 cm snowpacks at Barrow, Alaska (Tieszen 1974). Yet, an examination of snowpack observations from Daring Lake (Nobrega and Grogan 2007, Rees et al 2010 and surrounding regions (Derksen et al 2009) in relation to outputted thresholds from a penetration depth model by Zhou et al (2003) indicates that a portion of incoming light most likely penetrates thin (<5 cm) or patchy snowpacks during snowmelt at DL. Presently, subnivean photosynthesis has not been quantified over the vast majority of Arctic vegetation species or low Arctic sites. The explicit representation of subnivean photosynthesis within VPRM is therefore beyond the scope of this study, especially since GEE 0 describes snow and growing season GEE at both DL and IV with low RMSEs and MAEs.
Uncertainty in VPRM estimates of NEE was reduced when snow season respiration was calculated as a function of T soil (RESP s ) (table 2) Modeling subnivean respiration as a function of T soil rather than T air prevented the magnitude of respiration from being overestimated at the end of the snow season, when T air was consistently warmer than T soil (figures 3 and 5). This is consistent with previous findings that freeze-thaw temperature fluctuations accompanying snow melt do not substantially influence effluxes of CO 2 (Grogan et al 2004. Calculating subnivean respiration as a function of T soil allowed VPRM to simulate both the gradual, steady increase in respiration accompanying snowmelt, and midwinter declines in soil respiration (Bokhorst et al 2010).
However, when respiration was modeled year-round as a function of T soil (RESP Tsoil ), the RMSE and MAE errors were substantially larger than the errors in either RESP 0 or RESP s (table 3). The errors in RESP Tsoil were largest when SCA was <50%. As snow cover diminished, the contributions of aboveground biomass to respiration increased accordingly. The rate of respiration by vegetation is driven primarily by T air rather than T soil , and is therefore best modeled as a function of T air . As a result, unreliable estimates of growing season NEE would likely be generated using RESP Tsoil . Therefore, although estimates of subnivean respiration were improved when subnivean respiration was calculated as a function of T soil , uncertainty in VPRM estimates of NEE at the start and end of the growing seasons were reduced by calculating respiration as a function of T air (figure 6). These reductions in model uncertainty are important because the cold season is so long in the Arctic that even small biases in daily Table 2. Uncertainty in 1 May-30 June (MJ) and snow season (SS) estimates of NEE by VPRM both with ( s ) and without ( 0 ) representations of the influences of snow on respiration and GEE. Results are indicated for the Daring Lake MT calibration site, as well as the Ivotuk validation site, for years [2004][2005][2006][2007]. Mean absolute error (MAE) values are indicated first, followed by root mean squared error (RMSE) values in brackets. All error rates were calculated by comparing daily average eddy covariance NEE to daily average model NEE.  figure 7). Respiration was best simulated according to T air during the growing season because diurnal variability in T air influences the magnitude and timing of both photosynthesis and respiration. Conversely, the low thermal conductivity of an overlying snowpack decouples T soil and T air , and respiration persists throughout the snow season according to T soil . Due to the rapidity of this transition, and the fact that spring snowmelt usually occurs only once per year in Arctic regions, the use of a 50% threshold did not introduce any discontinuities in estimates of NEE. Differentiating the snow and growing seasons according to SCA, and simulating respiration as a function of the dominant temperature for each season, therefore allowed for reduced uncertainty in estimates of NEE at both sites (2004)(2005)(2006)(2007).

Conclusions
In Arctic regions, the timing of snow onset and melt influence the rates of photosynthesis and respiration. The importance of snow cover transitions for NEE suggests that insights into the northern carbon cycle and its response to changing snow conditions may be gained by representing the influence of snow on NEE. The feasibility of incorporating remote sensing observations of snow into models of NEE was demonstrated by findings showing: (1) good agreement between time-lapse camera (<10 m) and remote sensing estimates of SCA from  MT (left) and Ivotuk (right) as observed using the eddy covariance technique (black), and as estimated by the RESP 0 & GEE 0 (orange) and RESP s & GEE 0 (blue). Within each plot, the date where estimates from the two models appear to merge represents the day at which SCA initially decreases below 50%.
Landsat (30 m) and MODIS (500 m); and (2) associations between in situ NEE and SCA at Daring Lake, NWT (May-June 2010). Uncertainty in VPRM estimates of NEE at two low Arctic sites was reduced by representing the decoupling effects of a snowpack on T soil and T air . Estimating subnivean respiration as a function of T soil prevented respiration from being overestimated when it was limited by cool T soil at the start/end of the snow season, and enabled variability in cold season NEE to be simulated. The timing and magnitude of photosynthesis at the start and end of the snow season were best captured by GEE 0 , which used an implicit approach to simulate the influences of cold temperature, senescent vegetation and diminished sunlight on hindering photosynthesis. The resulting VPRM formulation, containing an implicit representation of the effects of SCA on photosynthesis and an explicit representation of the influence of SCA on respiration, had diminished RMSEs and MAEs across both sites and all years.
Climate change is predicted to both increase Arctic snow accumulation and diminish the length of the Arctic snow season (AMAP 2011). Previous studies at Daring Lake have found that natural inter-annual variability in snow melt timing did not markedly affect early or total growing season CO 2 flux (Humphreys and Lafleur 2011), but that plots with artificially increased snow depth and duration showed altered CO 2 fluxes upon snow melt (Buckeridge and Grogan 2010). Incorporating satellite observations of SCA into biospheric carbon flux models could therefore allow the snow and growing seasons to be delineated, enable snow season influences on respiration to be represented, and permit reduced uncertainty in estimates of the Arctic  carbon cycle. Insights could therefore be gained into the regional scale response of the Arctic carbon cycle to altered biological, meteorological and cryospheric conditions. Future work will consist of applying the snow season formulation to a variety of sub-Arctic, low Arctic and high Arctic sites in order to determine if model error may be reduced in all cases.