Decoupling between Plant Productivity and Growing Season Length under a Warming Climate in Canada’s Arctic ()
Received 9 May 2016; accepted 13 August 2016; published 16 August 2016
1. Introduction
Net primary productivity (NPP, g biomass m−2 y−1) is fundamental to an arctic tundra ecosystem’s functions and services, such as food for wildlife and feedback to climate [1] - [3] . Arctic temperatures have increased at almost twice the global average rate in the past 100 years and are predicted to continue at a higher rate than the rest of the planet [4] [5] . Along with rising temperature, an increase in growing season length (GSL, d) during the last few decades has also been reported [5] . Measurements of and GSL showed coincident declines from low latitudes to high northern latitudes [6] . Tundra warming experiments have also shown increases in both plant productivity and GSL [7] . Consequently, many climate-vegetation models have assumed a uniform response of plant growth to temperature and GSL, and reduction of arctic tundra [8] [9] .
A circumarctic study of relationships between shrub growth and climate variables (including GSL), however, raised doubt about this assumption [10] . Shrub growth at many Arctic sites, especially in the Canadian Arctic, was found to have no significant relationship with climatic variables at the 90% confidence level. In other words, shrub growth and GLS at these sites were decoupled.
Building on the findings of the paper [10] , we aim to investigate if such decoupling has also occurred for other tundra classes. Also, we will investigate the reasons of the decoupling. Among many possible reasons, we noticed that the distance between a shrub growth study site and the nearest climate station can be quite large in some cases, especially in the Canadian Arctic where climate stations were very sparse [10] . The difference between the climate conditions at a shrub site and those at its nearest climate station could thus weaken and even decouple the relationship. To circumvent this difficulty, in this study we use leaf NPP (NPPleaf) and GSL derived from the same source of long-term remote sensing time series over various tundra ecosystems in Canada’s Arctic.
2. Material
2.1. Study Areas
According to the Canadian ecozone classification [11] , two study areas (i.e., Ivvavik National Park (NP), Bathurst caribou summer range and calving ground) are located in the Southern Arctic Ecozone (Figure 1). The other two (i.e., Sirmilik and Torngat Mountains NP) in the Northern Arctic and Arctic Cordillera Ecozones.
The 1760 m high British Mountains dominate Ivvavik NP in northern Yukon. The park escaped glaciations
![]()
Figure 1. Location of four study areas in the Canadian Arctic, including Ivvavik, Sirmilik, Torngat Mountains National Parks, as well as the Bathurst caribou summer range and calving ground.
during the last Ice Age, except along its cottongrass tussock dominated northern coast. The Bathurst caribou habitat is a relatively flat rolling tundra shield, located in NWT and Nunavut. The Sirmilik NP is composed of the mountain and upland surrounding Oliver Sound, the rugged plateau of eastern Borden Peninsula, and Bylot Island. The Bylot Island features glaciers and a major seabird colony on its southwestern coast. Although its latitude is further south comparing to the Southern Arctic study areas, the Torngat Mountains NP also belongs to Northern Arctic and Arctic Cordillera Ecozones. The mean annual temperature and total precipitation at the Qikiqtarjuaq climate station, the nearest station from the Torngat Mountains NP in the same ecozone, indicate that it is colder and drier than the two Southern Arctic study areas (Table 1). The Torngat Mountains NP composes the George Plateau and the Torngat Mountains, which is the highest peak (1652 m) in east coast North America and dotted with remnant glaciers.
2.2. Data Sources
To produce NPPleaf and GLS for each tundra class in the four study areas from 1985 to 2013, we used long-term satellite data the Advanced Very-High-Resolution Radiometer (AVHRR). The AVHRR data are pre-processed 10-day composites at 1-km spatial resolution, including the red and near-infrared surface reflectance (i.e., ρr and ρnir) and cloud probability. Pre-press procedures include geo-referencing, the Bidirectional Reflectance Distribution Function (BRDF) correction, atmospheric correction, cloud indexing, and inter-sensor normalization [12] .
Field measurements for calibrating and validating the remote sensing products of NPPleaf and GLS include leaf biomass and percentage vegetation cover. Leaf biomass and percentage vegetation cover were collected during July and August at 11 sites in the Ivvavik NP in 2008, at 16 sites in Torngat Mountains NP in 2008, at 11 sites in Sirmilik NP in 2010, and at 34 sites in the Bathurst caribou habitat during 2005, 2013 and 2014.
To directly match an AVHRR pixel, a suitable site would have to be homogeneous and at least 3 km by 3 km because of the possibilities of different orientations and partial overlapping between them. Our experience suggested that it could be difficult to find enough such sites in the Arctic, where topographical variations often resulted in a heterogeneous landscape. Therefore, our criteria for site selection were at least 90 m by 90 m or the size of three Landsat pixels. The 30-m resolution Landsat images were then used to bridge the scale difference between leaf biomass measurement sites and AVHRR pixels. We developed a Landsat mosaic for each of the four study areas cloud-free circa 2000 scenes, following paper [13] . Given the mismatch in the dates often found between field measurements and their corresponding Landsat mosaic, a correction is needed to convert the vegetation indices calculated using the Landsat mosaic to that on the dates of field measurements. To carry out the correction, we used the daily 250-m data of the Moderate Resolution Imaging Spectroradiometer (MODIS), which acquires data daily approximately 15 min after Landsat data in the same polar orbit.
Tundra classes at 1-km resolution for each of the four study area were aggregated from 30-m Landsat-based land cover maps [14] - [16] . Soil moisture regime index was determined using field survey records and DEM data [17] . The DEM tiles at 1:50,000 were downloaded from Geobase (www.geobase.ca). The climate data (temperature and precipitation) are available from Canadian Daily Climate Data (CDCD) (www.climate.weatheroffice.gc.ca).
3. Methods
3.1. In-Situ Measurements of Vegetation Variables
Growing season lengths in the Canadian Arctic are very short, ranging from a few weeks to a few months. Dur-
![]()
Table 1. Study areas, their attributes, and the nearest climate station with the good data record in the same ecozone.
ing the short growing season, leaf biomass values are relatively stable from middle July to middle August. This period was thus selected as the most suitable period for conducting field measurement of leaf biomass, for the purpose of calibrating and validating remote sensing products.
We select a field measurement site to be relatively homogenous and of a minimum size of 90-m × 90-m, to match with a Landsat pixel [18] . At each site, we typically sampled five to twenty 1-m × 1-m plots [18] . The number of plots per site was determined by balancing the need to go to more sites in a study area with the need to have more plots in a site. Statistics from our results indicate that to ensure a site is well represented (i.e., the sampling error <20%), we need to have five plots per site for sites with vascular plant percentage cover >20%, and 20 plots per site for sites with lower vascular plant percentage covers. Field measurements in the remote Canadian Arctic are expensive and time-consuming because many sites are accessible only by helicopter. To maximize the value of helicopter time, we usually included measurements other than leaf biomass and percentage covers (e.g., permafrost depth, soil properties and moisture regime, roots biomass, and foliage photosynthesis capacity) in a field campaign. Consequently, a team of 4 to 5 researchers usually was only able to measure leaf biomass and percentage covers for one or two sites per day. A summer’s field campaign thus usually allowed us to collect leaf biomass and percentage covers 10 - 25 sites in a study area. To avoid sites overly concentrated at the dominant land cover classes, we intentionally included sites with at the high and low ends of biomass range.
At each plot, percentage covers of vascular plant species were visually estimated in the field and corrected by digital photo classification in the laboratory [19] . All plants were then harvested, identified to species, sorted into dead versus live as well as leaves versus stems, and weighed in the field. Samples of these leaves and stems were also taken to the laboratory and oven-dried and weighed to obtain the oven-dry leaf biomass. We calculated the values of leaf biomass and percentage cover at each site as the average of all plots at the site.
3.2. Method for Quantifying NPPleaf
NPP includes the production of leaves, stems, roots, flowers, and fruits of vascular plants, as well as the production of non-vascular plants (e.g. lichen, mosses). For ecosystems in the Canadian Arctic, it was suggested that belowground NPP may strongly related to its aboveground counterpart [20] , although paper [21] found there were not enough root production measurement to be conclusive. Therefore, NPPleaf can be used as a good surrogate for NPP in Canada’s Arctic.
By definition, NPPleaf is the net increase of leaf biomass from the start of growing season to its peak value. As such, NPPleaf is given by (the seasonal peak SRVI ? the SRVI at the SOS) × the slope of the relationship between leaf biomass and SRVI. Here SRVI is the AVHRR-based simple ratio vegetation index (=ρnir//ρr). In this study, we used SRVI, instead of the commonly used normalized differential vegetation index (NDVI), because the relationships between NDVI and leaf biomass over arctic tundra ecosystems are usually nonlinear and have a lower explanation power than that between SRVI and leaf biomass [22] .
Despite the selection of the best day within a 10-d period through compositing techniques and other pre- processing efforts, the remnant cloud contamination and aerosol variations could still result in an up to 40% error in the seasonal peak value of SRVI [13] . To reduce the large error in the seasonal peak value of SRVI, we constructed seasonal profiles of SRVI for each land cover class, instead of each AVHRR pixel [23] . We divided the pixels of a tundra class in a given 10-d composite period into four categories: clear sky, lightly, moderately, and heavily cloud contaminated, with the cloudiness index < 20%, 20 to <40%, 40 to <60%, and ≥ 60%, respectively. SRVI over cloud contaminated pixels in a tundra class then was corrected against that of clear-sky pixels in the class during the same 10-d composite. Table 2 lists the multipliers for correlating (the value of mean SRVI-1) of cloud contaminated pixels in a tundra class to that of clear pixels in the class for the four study areas in the Canadian Arctic. With these correction coefficients, we were able to reduce the error in SRVI of a land class to <10%.
Next, we determined the slope of the relationship between leaf biomass and SRVI while describing the quantification of the SRVI at the SOS in section 3.3. Due to the spatial mismatch between the size of a field measurement site and that of an AVHRR pixel, we cannot directly develop the relationship between leaf biomass and AVHRR-derived SRVI. Instead, we used Landsat mosaic to bridge their scale difference. Landsat images were pre-processed, which included 1) co-registration of images; 2) atmospheric correction to obtain surface reflectance in red and near infrared, using a procedure based on 6S (Second Simulation of Satellite Signal in the Solar Spectrum) described by reference [24] . The atmospheric parameters (i.e., atmospheric water vapor mean,
![]()
Table 2. Multipliers for correlating (the value of mean SRVI -1) of cloud contaminated pixels in a tundra class to that of clear pixels in the class for the four study areas in the Canadian Arctic.
total ozone mean, and optical depth and ocean Mean) were developed from the MODIS instrument, which acquires data daily approximately 15 min after Landsat data in the same polar orbit; and 3) cloud and shadow masking. For creating the mosaic, we kept a baseline image intact and normalized other images to it through equations developed over the overlapping areas.
Because of the 16-day repeat cycle and frequent cloud cover, acquiring Landsat images exactly corresponding to the field measurement date is often impossible. As a result, substitute Landsat images acquired on another date or year are often used. Given that there are substantial seasonal and inter-annual variations in foliage biomass, such substitution can significantly increase the error in the relationship between leaf biomass and Landsat-derived signals due to the mismatch between the date of Landsat image acquisition and that of field measurement. To minimize this error, we made a correction using MODIS daily 500-m surface reflectance data. We developed relationships for MODIS band 1 (or 2) between the date of Landsat acquisition and that of field measurement and then applied this relationship to the Landsat red (or near-infrared) surface reflectance to correct the Landsat reflectance to the day of field measurement, using an approach similar to paper [25] . Table 3 lists the statistics for the relationships between leaf biomass and Landsat-derived SRVI for the four study areas in the Canadian Arctic. Applying these relationships to the Landsat mosaic, we developed leaf biomass maps for each study area, as exemplified in Figure 2 for the Ivvavik NP.
![]()
Table 3. Statistics for the relationships between leaf biomass and SRVI derived from Landsat imagery or AVHRR data for the four study areas in the Canadian Arctic. Here n is the sample size.
![]()
Figure 2. Landsat-derived leaf biomass map over the Ivvavik National Park. We used DEM data for a three-dimensional effect.
Using these leaf biomass maps, we can calculate 1-km resolution leaf biomass distribution and relate them to AVHRR-derived SRVI. To reduce the random error in AVHRR SRVI for a clear sky pixel, we averaged AVHRR-derived SRVI and leaf biomass to a window of 25 pixels, as a compromise between final accuracy in the averaged AVHRR SRVI and a large enough sample number for the regression [13] . The slopes of the relationships between leaf biomass and AVHRR-derived SRVI for the four study areas in the Canadian Arctic was also listed in Table 3.
3.3. Method for Determining Plant Seasonality
We used the Biophysically-based and Objective Satellite Seasonality Observation Method (BLOSSOM) for determining GSL [26] . Biophysically, the date of SOS or EOS can be defined as the day of the year on which the leaf biomass of deciduous shrubs and herbs equals zero [26] . Due to the existence of evergreen shrubs in most tundra sites, AVHRR-derived SRVI threshold at the SOS (or EOS) was determined by considering contributions from both evergreen shrubs and other components. The other land surface components might include dead leaves and stems of deciduous shrubs and herbs, lichens, mosses, bare soil, rocks, and shadows. The contribution of evergreen shrubs to the SOS threshold was determined by the percentage cover of evergreen shrubs in a tundra class and their surface reflectance while the contribution from other components was given by (-intercept/ slope) for each study area listed in Table 3. Table 4 lists the percentage cover of evergreen shrubs in a tundra class and the SOS (or EOS) AVHRR-derived SRVI threshold for each class in the four study areas.
As shown in Figure 3, the SOS was determined as the day of the year on which class average AVHRR SRVI first became larger than the threshold in the spring. Similarly, the EOS was determined as the day of the year on which class average AVHRR SRVI first became less than the threshold in the fall. Since snow could sometimes revisit a vegetated area around SOS or EOS in the Arctic, and result in multiple crossings of AVHRR SRVI through the threshold line, only the first crossing with certainty in the spring would be determined as SOS. All other days before this date should have no class mean AVHRR SRVI―its one standard estimation error were larger the threshold. Similarly, only the first crossing with certainty in the fall would be determined as EOS. The GSL was then determined as the difference between the date of EOS and the date of SOS plus one day.
3.4. Calculation of Soil Moisture Regime Index and Coupling Strength of a Relationship
Paper [17] determined soil moisture regime index (SMRI) for the Firth River Valley, a part of the Ivvavik NP, by information of soil texture and topography. They set the value of SMRI for a given soil moisture regime to be: 1 stands for the driest xeric regime, 2 for subxeric, 3 for submesic, 4 for mesic, 5 for subhygric, 6 for hygric, 7 for subhydric, and 8 for the wettest hydric regime. Following the same method, we calculated the SMRI for each class in the four study areas, by combining information from soil texture survey during our field campaigns and
![]()
Figure 3. An example showing the determination of SOS, EOS, and GSL for the barren class in Sirmilik NP in 2007. The vertical bars show one standard error in AVHRR SRVI of the class while the horizontal bars show the corresponding under (or over-) estimation of SOS and EOS.
![]()
Table 4. Percentage cover of evergreen shrubs and SOS (or EOS) threshold for each class in the four study areas. Also listed are the mean values of NPPleaf and GSL, long-term trends of NPPleaf and GSL during 1985-2013, the correlation coefficient Rcs between NPPleaf and GSL, Rdt between de-trended NPPleaf and de-trended GSL, soil moisture regime index, and the percentage of a tundra class out of the total tundra-covered area in a study area. Long-term trends statistically significant at the 99%, 95%, and 90% levels are labeled as ***, **, and *, respectively.
topographic information from DEM data (Table 4).
We defined the coupling strength between NPPleaf and GLS as their correlation coefficient (Rcs). To better understand the causes of the decoupling, we further decomposed the relationship into two components: the correspondence of inter-annual variations and the agreement of long-term trends. The pure inter-annual variation of NPPleaf (or GLS) refers to the change from year to year after the long-term trend has been removed. Therefore, we calculated the contribution of the correspondence between interannual variations to the coupling strength between NPPleaf and GLS as the correlation coefficient (Rdt) between the de-trended NPPleaf and de-trended GLS. The de-trended NPPleaf (or GLS) was calculated as NPPleaf (or GSL) ? its trend value at a given year. Since there is only one value for the long-term trend of NPPleaf and another for that of GLS, direct calculation of the contribution of the agreement between long-term trends to the coupling strength between NPPleaf and GLS is not possible. Instead, we can estimate the contribution of the agreement between long-term trends as (Rcs − Rdt).
4. Results and Discussions
4.1. States and Trends of NPPleaf and GSL across the Four Study Areas
The mean NPPleaf of a tundra class during 1985-2013 ranged from 23.5 to 80.3 g∙m−2∙y−1 over the two southern study areas, and from 4.9 to 72.6 g∙m−2∙y−1 over the two southern study areas (Table 4). Our results are comparable with results of other studies in the Canadian Arctic, with aboveground NPP (ANPP) ranging from 0.7 to 60 g∙m−2∙y−1 for high Arctic sites and 28 to 125 g∙m−2∙y−1 for low Arctic sites [27] . Note that NPPleaf is the largest component of ANPP, we expect NPPleaf being somewhat lower than ANPP.
Within a study area, the values of mean NPPleaf of a tundra class during 1985-2013 could differ by several folds. For example in the Ivvavik NP, the willow-sedge pediment drainage channel class (i.e., I20) had a mean NPPleaf during 1985-2013 of 70 g∙m−2∙y−1, contrasting to that of the hedysar-avens inactive alluvial terrace class (I23) at 33 g∙m−2∙y−1. The largest difference in NPPleaf was found between the shrub-herb-lichen-bare class (T24) and the rock outcrop low vegetation cover class (T38) in the Torngat Mountains NP: 72.8 g∙m−2∙y−1 versus 4.9 g∙m−2∙y−1. In comparison, the difference in GSL among classes within a study area was not as large as that in NPPleaf. For example, the GSL of the willow-sedge pediment drainage channel class was 130 d, while that of the hedysar-avens inactive alluvial terrace class 110 d in the Ivvavik NP. The largest difference in GSL was found between the herb-shrub class (T23) and the rock outcrop low vegetation cover class (T38) in the Torngat Mountains NP: 127 d versus 52 d.
When averaged over all tundra classes in a study area weighted by area percentages, the mean values of NPPleaf of the two southern Arctic study areas (i.e., Ivvavik NP and the Bathurst caribou summer range and calving ground) during 1985-2013 were, respectively, 55.5 g∙m−2∙y−1 and 56.5 g∙m−2∙y−1. For the two northern Arctic study areas, the mean values of NPPleaf during the same period were much smaller: 13.4 g∙m−2∙y−1 over the Sirmilik NP and 14.4 g∙m−2∙y−1 over the Torngat Mountains NP. Correspondingly, the values of study area averaged GSL during 1985-2013 in the southern Arctic study areas were 127 d over the Ivvavik NP and 123 d over the Bathurst caribou habitat. To the contrary, the values of study area averaged GSL in the northern Arctic study areas were 71 d over the Sirmilik NP and 81 d over the Torngat NP.
Note that the percentage of a tundra class was calculated against the total tundra-covered area in a study area. Glaciers, water bodies, and wetland classes were excluded in the calculation. In the Sirmilik NP, the glaciers account for about 40% of total area of park land. If they were included in the study area averaged NPPleaf and GSL calculation, their values would be even smaller. We also noticed that these study area averaged values of NPPleaf and GSL were mainly controlled by their most dominant classes. For example, in the Torngat Mountains NP, the low vegetation cover (bare soil, rock outcrop) class alone accounted for 66.4% of total tundra area in the park, followed by the lichen barren at 19.3%. Small classes, such as the deciduous shrub class (0.4%) and the shrub-herb-lichen-bare class (0.2%), did not contribute to the study area averaged NPPleaf despite their much higher values of NPPleaf (Table 4).
As for the long-term trends of NPPleaf and GSL, large variations were also observed among tundra classes within a study area. All classes in the four study areas were found to have increasing trends in their NPPleaf and GSL, although some of the trends were not statistically significant at the 90% confidence level. The largest difference in the long-term trends of NPPleaf was observed between the rock outcrop low vegetation cover class (T38) and the shrub-herb-lichen-bare (T24): 0.29 versus 2.74 g∙m−2∙y−1 per year. Typically, the difference in the long-term trends of GSL between classes in a study area was smaller than that of NPPleaf, which was consistent with the finding for their mean values. More analyses of these long-term trends are reported in section 4.4.
4.2. Relationships between NPPleaf and GSL during 1985-2013
Figure 4 shows examples of covariations between NPPleaf and GSL. The NPPleaf of the willow-birch moist slope class in the Ivvavik NP increased significantly during 1985-2013, with R2 = 0.25, p = 0.005, and sample size n = 29. During the same period, its GSL also increased significantly (R2 = 0.15, p = 0.04). Yet, despite of their common increasing trends, they were decoupled (i.e., not significant at the 90% confidence level), with R2 = 0.06, p = 0.21, and n = 29.
Seven out of the 11 tundra classes (i.e., 64%) in the Ivvavik NP, their NPPleaf and GSL were decoupled during 1985-2013 (Table 4). Similar to the finding of paper [10] , many of shrub classes in the Ivvavik NP (e.g., class I2 the willow-horsetail wet slope, I4 the willow-birch moist slope, I20 the Willow-sedge pediment drainage channel, and I25 the willow-coltsfoot drainage channel) had decoupled relationship between NPPleaf and GSL during 1985-2013. Other shrub classes (e.g., class I10 the willow floodplain) had a strongly coupled relationship between NPPleaf and GSL during the same period. The results were similar for other land cover classes. For example, classes that had decoupled relationship between NPPleaf and GSL include I1 the alpine slope, I18 the cottongrass tussock, and I26 the sedge tussock. On the other hand, other classes (e.g., I3 the rock lichen, I22 the sand/silt, and I23 the hedysar-avens inactive alluvial terrace) had a coupled relationship between NPPleaf and GSL, although their coupling strengths were not as significant.
For the Bathurst caribou habitat, similar variable results were found. Four out of 9 tundra classes (i.e., 44%) in the Bathurst caribou habitat, their NPPleaf and GSL were decoupled during 1985-2013 (Table 4). Both shrub classes (i.e., B16 and B19) had a decoupled relationship between NPPleaf and GSL. Different from the Ivvavi NP, the herb class in the Bathurst caribou habitat had a coupled relationship between NPPleaf and GSL during 1985-2013. As for the northern arctic Sirmilik study area, NPPleaf and GSL were strongly coupled for all tundra classes (Table 4), including the barren class (Figure 4).
To further investigate the impacts of climate change on the relationship between NPPleaf and GSL, we compared their coupling strength between the first 10-year period (i.e., 1985-1994) and the last 10-year period (i.e., 2004-2013). For the two southern Canadian Arctic study areas, we found the majority of classes (75% for the Ivvavik NP and 78% for the Bathurst study area, respectively) had a reduced coupling strength from the first 10-year period to the last 10-year period. To the contrary for the two northern Canadian Arctic study areas, we found only about half of classes (67% for the Similik NP and 38% for the Torngat study area, respectively) had a reduced coupling strength between the two periods.
A question then arises as to what controlled the coupling strength of the relationship between NPPleaf and GSL? We addressed this question by investigating both components of the relationship: the correspondence of inter-annual variations and the agreement of long-term trends, in the flowing two sections.
![]()
Figure 4. Trends in NPPleaf and GSL of the willow birch moist slope class in Ivvavik NP, and the barren class in Sirmilik NP during 1985-2013. Error bars show uncertainties at one standard deviation.
4.3. Correspondences between Interannual Variations in NPPleaf and GSL
We analyzed the correspondence between their inter-annual variations against various mean characteristics of a land cover class within a study area, such as long-term mean leaf biomass and GSL during 1985-2013, mean elevation of all pixels of the class, and SMRI. We found that the mean GSL to be the most significant variable for explaining the change in the correspondence.
As Figure 5 shows, the ratio of sign agreement between de-trended NPPleaf (i.e., NPPleaf ? trend value of NPPleaf) and de-trended GSL (i.e., GSL ? trend value of GSL) for the barren class in the Sirmilik NP during 1985-2013 was 100% for years with GSL < 30 d. The ratio of sign agreement was computed as years with agreed signs divided by the total years in a given GSL range. The high ratio values for years with GSL < 30 d are understandable because it usually takes a couple of weeks to dry out soil, even for sites with very coarse soil texture in the Arctic [20] . Because leaf biomass usually peaks near the middle of a growing season, there would be almost no chance of drought for years with GSL < 30 d. The ratio decreased quickly for years with GSL > 30 d to ~50%, the threshold corresponding to complete decoupling. In other words, in years with longer GSL, factors such as soil moisture limitation, herbivore grazing, or insect defoliation had a much bigger chance to become influential, resulting in decoupled NPPleaf and GSL. The ratio of sign agreement in years with GSL > 30 d in Ivvavik NP was slightly higher than that in Sirmilik NP, probably because Ivvavik NP has higher precipitation and thus might take a longer period to develop drought.
The overall lower ratio of sign agreement for the willow-birch moist slope class in the Ivvavik NP might also be explained by the difference in temperature sensitivity for NPPleaf and GSL. During 1985-2013, GSL of the class was significantly correlated with the mean temperature in the months of May, June, September, and October (i.e., spring and fall temperature in the southern Canadian Arctic), with R2 = 0.31 and p < 0.01. However, they were not correlated with mean summer temperature during July and August, with R2 = 0.006 and p > 0.1. On the other hand, NPPleaf was related to the July and August mean temperature (R2 = 0.12 and p < 0.1), but not with spring and fall temperature (R2 = 0.0002 and p > 0.1), similar to the finding of by paper [28] . No significant relationship was found between precipitation and NPPleaf or GSL for all time periods.
These different sensitivities, in fact, support the argument that other factors such as soil moisture limitation, herbivore grazing, or insect defoliation might constrain plant growth. Otherwise, if all these factors’ impacts were negligible, then an earlier and better spring growth should have also resulted in higher NPPleaf.
For the barren class in the Sirmilik NP, NPPleaf was decoupled with the mean July and August temperature (R2 = 0.093, p = 0.107, and n = 29). The decoupled relationship was a remarkable contrast to the high coupling
![]()
Figure 5. Percentage of sign agreement between de-trended NPPleaf and de-trended GSL in a GSL range for each class in the Sirmilik NP during 1985-2013. See the class label in Table 2.
strength between its NPPleaf and GSL (R2 = 0.34 and p < 0.01) for the class. We suspected that the tempo-spatial mismatch between climate and NPP data might be responsible for the contradiction.
The percentage of years with GSL < 30 d increases as the mean GSL of a tundra class decreases. Therefore, we used the mean GSL of the class as a proxy to represent the chance for factors such as soil moisture limitation, herbivore grazing, and insect defoliation to become influential. A significant relationship was found between the coupling strength between de-trended NPPleaf and de-trended GSL (Rdt) and the mean GSL of a tundra class (
) over the four study areas (Figure 6):
(1)
As shown in Figure 6, the inter-annual correspondence between de-trended NPPleaf and de-trended GSL became statistically insignificant for some of the classes that had a long-term mean GSL > 100 d during 1985-2013. When the long-term mean GSL was > 120 d, nearly all classes had an insignificant correspondence between de-trended NPPleaf and de-trended GSL.
4.4. Agreements between Long-Term Trends in NPPleaf and GSL
As for the agreement between long-term trends, we found that the SMRI was the dominant factor. Plotting the SMRI against the contribution of the agreement between long-term trends of NPPleaf and GSL to the coupling strength (i.e., Rcs ? Rdt), we found a significant reduction of this contribution when moving from the wettest tundra class to the driest tundra class for all four study areas (Figure 7, Table 5). For some of the driest tundra ecosystems, such as the rock outcrop low vegetation cover in the Torngat Mountains NP or the low vegetation cover (bare soil, rock outcrop) in the Bathurst caribou summer range and calving ground, this contribution even became negative.
These results clearly demonstrated the limitation of soil moisture on the contribution of NPPleaf and GSL trends to the coupling strength. One possible explanation is that a dry class with lower SMRI tended to limit the increase in NPPleaf (Figure 8). The ratios of NPPleaf trend of the wettest class to that of driest class during 1985-2013 were 4.8, 7.0, 2.7, and 6.2 for Ivvavik, Bathurst, Sirmilik, and Torngat study areas, respectively. To the contrary, GSL trends tended to increase for dry classes, probably because of the effect of enhanced temperature increase at higher and usually drier locations [29] . The higher GSL increase rate together with lower NPPleaf increase trend at the dry classes resulted in poor correspondence between NPPleaf and GSL, and thus had a minimal or even negative contribution to their coupling strength.
Combining Equation (1) and that in Table 5 for the influence of SMRI of the class, we estimated the coupling strength for each class in the four study areas:
![]()
Figure 6. The relationship between the mean GSL of a tundra class and the correlation coefficient between de-trended NPPleaf and de-trended GSL (Rdt) over the four study areas during 1985-2013.
![]()
Figure 7. SMRI against the contribution of agreement between long-term trends of NPPleaf and GSL to the coupling strength (i.e., Rcs ? Rdt) for each of the four study areas. Table 5 lists statistics for the relationships.
![]()
Table 5. Statistics for the linear relationships between the SMRI and (1) the contribution of the agreement between long-term trends of NPPleaf and GSL to their coupling strength, as calculated by Rcs - Rdt; (2) the GSL trend; or (3) the NPPleaf trend for each study area during 1985-2013.
![]()
Figure 8. The relationship between SMRI and trends of NPPleaf and GSL during 1985-2013 for each of the four study areas. Table 5 lists statistics for the relationships.
![]()
(2)
with R2 = 0.7, p < 0.01, and n = 36. Together, the ![]()
and SMRI can explain 70% of the variation in the coupling strength among all tundra classes in the four Canadian Arctic study areas (Figure 9). According to Equation (2), the relationship between NPPleaf and GSL became decoupled when ![]()
was >136 d for tundra classes of SMRI = 6.5, but only when ![]()
was >118 d for tundra classes of SMRI = 1 in the Bathurst caribou summr range and calving ground. In comparison for the Sirmilik NP, the relationship between NPPleaf and GSL became decoupled when ![]()
was >123 d (or 116 d) for tundra classes of SMRI = 6.5 (or 1). Similar results were also found for the Ivvavik NP and Torngat Mountains NP. These results thus predict that under a warming climate the decoupling could occur when the ![]()
would have increased to >116 d for dry tundra classes and >123 d for moist tundra classes in the northern Arctic. In the southern Arctic, the decoupling could occur when the ![]()
would have increased to >118 d for dry tundra classes and >136 d for moist tundra classes.
5. Conclusions
We estimated NPPleaf and GSL during 1985-2013 for all tundra classes in two northern Canadian Arctic study areas and two Canadian southern Arctic study areas, using long-term remote sensing data and field vegetation measurements.
Our results showed that 44% - 64% of tundra classes in the two southern Arctic study areas had a decoupled NPPleaf and GSL relationship during 1985-2013. The decoupling occurred for both shrub classes and non-shrub tundra classes. In the two northern Arctic and Arctic Cordillera study areas, we didn’t find decoupling for all tundra classes.
The strength of coupling increased significantly with the reduction of mean GSL of a tundra class, probably because a very short growing season giving less chance for other factors (e.g. soil moisture limitation) to become influential. The mean GSL of a tundra class thus might be a good indicator of the chance that these factors
![]()
Figure 9. Comparison of coupling strength Rcs between NPPleaf and GSL estimated using the long-term mean GSL and SMRI against that measured directly from NPPleaf and GSL for each class during in the four study areas.
will become influential. However, a further test of this hypothesis is needed, and collection of these datasets, along with NPPleaf and GSL is thus of great importance.
The concurred increasing trends of NPPleaf and GSL enhanced the coupling strength for most tundra classes. The enhancement, however, decreased significantly for drier classes, for which the long-term trend of NPPleaf was often insignificant due to soil moisture limitation, whereas GSL increased at a similar or even higher rate than wetter classes.
Together, the mean GSL and SMRI of a tundra class can explain 70% of the variation in the coupling strength among all tundra classes in the four Canadian Arctic study areas during 1985-2013. This result suggests that the decoupling between plant productivity and GSL would be more likely to occur under a warming climate if the tundra class had a mean GSL > 116 (or 123) days with a dry (moist) soil moisture regime.
Acknowledgements
The study was funded by the Canadian Space Agency’s Government Related Initiatives Program (GRIP), and the NWT Cumulative Impact Monitoring Program (CIMP). We thank the contribution of many northerners to the field measurements, including Alexander Gordon, Jayneta Pascal, and Kayla Arey from Aklavik, NWT; Simeonie and Jonas Johncena of Pond Inlet, Nunavut; Roy Judas and Brain Kodzin from Wekweeti. We appreciate the constructive comments by anonymous reviewers.