2.1 Climate Change

Decadal changes in temperature of the near-surface atmosphere in the Arctic have profound implications for the loss of snow and ice and thus their feedback to regional and global climate (Hinzman et al., 2013). Surface temperature changes are related to the air energy content (Graversen et al., 2008), which is an approximation of air heat content (i.e., the sum of enthalpy and latent heat, which are the respective functions of air temperature and humidity) (Pielke et al., 2004). Due to profound effects of the latter on the various dynamics, we use heat content as the metric of “environmental conditions.” We consider this climactic driver a primary forcing factor of change in the abiotic, biotic, and human systems.

An increase in near-surface heat can be conveyed via weather, event-scale impacts due to heat waves characterized by peak temperature and humidity as well as duration above a threshold. Gradual increase in the warming and duration of the above-freezing period leads to climate-scale changes, such as the timing of season onset and termination, their average heat content and duration. We consider changes for both types as external, that is, without accounting for how the Arctic land-surface will feedback to them. Even with this simplified view, we can distinguish two concepts. First (a), there can be temporal persistence of processes triggered by both pulse-scale and climate-scale changes (e.g., a brief heat wave may trigger an over-a-threshold behavior with long-term consequences; likewise, gradual changes in annual seasonalities may cause incremental but continuous changes in the Arctic system under consideration). Second (b), the impacts of heat content changes are often overlapping, and may lead to multiple feed-forward and feedback loops in the systems of impact, among which we target those that have longer-term implications. Arguably, mechanistic, process-level understanding of linkages generated in (a) and (b) constitute the core of Arctic climate change impacts that are of interest to science and society (Hinzman et al., 2013; Meier et al., 2014; Overland et al., 2016; Walsh et al., 2020). We consider the complexity of relevant pathways by using both temporal scales associated with an increase in the atmospheric heat content and consider questions that require a convergence science approach.

2.2 Industrialization

The second stressor about which we are concerned is industrial development. Globally, industrialization is defined as a period of social and economic change during which people's means of gaining subsistence shifts to minimize human drudgery and improve predictability of resource availability via systematization and simplification of processes, an extensive division of labor, substitution of mechanical for human energy, and replacement of small, localized, and uncertain sources of supply by large, networked and controllable ones (Shimkin, 1952). Industrialization is a complex phenomenon with many regional and temporal elements that result in particular histories and complex constellations of identities and socio-political groupings. Industrialization is also an accelerator of acculturation that has affected every society across the globe, as people in less industrialized societies borrow or adapt to features of more industrialized cultures. Industrial societies have all experienced dramatic increases in the per capita production of food, services and goods through the mechanization of manufacturing and agriculture. Industrialization depends on the social systems that provide labor (Robertson, 1991). A social system consists of individual human beings interacting with one another within certain continuing associations and institutions.

Russian industrialization of the Arctic is characterized by large-scale operations associated with the exploitation of non-renewable resources, and the construction of cities built around extraction cites and transportation networks (Zamyatina, 2023). This was accompanied by a large-scale influx of workforce who settled first temporarily, then permanently in such cities (Bolotova & Stammler, 2010). Alongside the development of non-renewable resources, indigenous peoples across Siberia, the Russian Arctic, and the Far East were incorporated into the state agricultural system, first through trading cooperatives linked to state procurement agencies, and then in the 1930s into collective farms (kolkhoz), and 1960s state farms (sovkhoz and gospromkhoz) for exploitation of renewable resources (mainly reindeer, fish, and fur) (Ziker, 2002). Oil and gas resource development, intensifying after the breakup of the USSR, is more distributed and reliant on the construction of linear infrastructures, such as pipelines, and exploitation of shift-worker regimes, rather than construction of industrial cities (Saxinger, 2016). Although recent expansion of oil/gas and mining associated infrastructure in the arctic has mostly occurred in Russia, there are also significant developments in the US and Canada (Bartsch et al., 2021).

3.1 Eco-Gradient

The Yamal Peninsula, West Siberia, Russia represents a microcosm of the changing Arctic, where the two novel stressors described in Section 2 interact with a dynamic social ecological system. The peninsula is a clearly bounded natural laboratory with a biogeographic gradient from the forest-tundra ecotone to the high Arctic (Figures S1a–S1c in Supporting Information S1) and extends over 700 km south-north (240 km east-west) from the northern terminus of the Polar Urals to the Kara sea, presenting four of the five Arctic bioclimatic subzones, from subzones E, D, and C in the main land of the peninsula, to subzone B in the adjacent Belyi Island (Walker et al., 2005). Yamal is predominantly underlain by the continuous and in the south by discontinuous permafrost (Figure S1b in Supporting Information S1). Yamal evinces variation in both plant and animal communities along a latitudinal gradient. For vegetation, there is a general decrease in productivity, height and biodiversity of plant types; and there is an increase in the ratio of mosses to vascular plants from south to north (Walker et al., 2009). There is also a decreasing number of animal species from south to north (Ryzhanovsky et al., 2016). These plant and animal gradients mirror those seen elsewhere in the Arctic across bioclimatic subzones.

3.2 Rapid Climate Change

The region is a hotspot of surface air temperature warming: June–July warming over the period 1991–2020 has led to an increase of +1.32°C as compared to the climate normal period of 1961–1990, or +2.02°C, relative to preindustrial levels (1850–1900), far exceeding the range of natural climate variability (Hantemirov et al., 2022). Mean annual temperature changes from 1961 to 1990 to 1991–2020 are about +1.5°C (Malkova et al., 2022). There is a positive trend in liquid precipitation (mean total annual is 390 mm, 2000–2019), accompanied by a decrease in snowfall (mean total is 223 mm), and an increased likelihood of rain-on-snow events (Forbes et al., 2016). Rapid warming increases the vulnerability of the permafrost to thawing (Malkova et al., 2022; Shpolyanskaya et al., 2022). Borehole measurements indicate one of the fastest warming rates of ground temperatures across the Arctic regions with continuous permafrost (Biskaborn et al., 2019), although this can vary with the type of landscape (Kaverin et al., 2017). As in much of the Arctic, warming has been linked to an increase in height and abundance of tall shrubs and a shift of the forest-tundra ecotone in the southern half of the peninsula (Frost & Epstein, 2014; Hantemirov et al., 2008; Mazepa, 2005; Shiyatov et al., 2007).

Considerably less data have been published to date concerning responses of terrestrial fauna in Yamal to climate change. However, studies do suggest that animal species' ranges have shifted northward with climate change, leading to the “borealization” of small rodent and bird (Sokolov et al., 2012) communities, and expansion of breeding ranges of red foxes (Vulpes vulpes L.) and hooded crows (Corvus cornix L.) (Sokolov, Sokolov, & Dixon, 2016). These have important consequences for food web structure and functioning (Ims et al., 2013a, 2013b).

3.3 The Built Environment

The built environment on Yamal includes towns and villages developed over the last 90 years, as well as a network of trading posts (faktoria), industrial extraction sites, and slaughterhouses. Linear infrastructures include a railway, and some short concrete roads. Long distance roads are maintained as ice-roads during the winter season. Industrial development has increased rapidly since the 1990s (Stammler, 2011). The largest industrial facilities situated in Yamal are Bovanenkovo, Sabetta, Noviy Port and Kharasavey, with several thousand workers each. The footprint of infrastructure and urban development in Yamal is not large in a spatial context, but much of the new infrastructure is linear and connects previously remote and relatively isolated communities (Kumpula et al., 2012). The 572 km Obskaya–Bovanenkovo railway (Figures S1b–S1c in Supporting Information S1), completed in 2011, is the world's northernmost railway (Terekhina & Volkovitskiy, 2020). Industrial activity brings large numbers of shift workers, commuting between urban centers in Yamal-Nenets Autonomous Okrug (YaNAO) and other cities of Russia and fly-in/fly-out settlements, such as Sabetta and Bovanenkovo. Industrial commuters, thus now outnumber the permanent population of native villages. Important for this paper are new stressors associated with industrialization, specifically the growing presence of industrial activities, such as construction in previously hard-to-access places, maintenance and transport of equipment and supplies for the gas industry, resulting in the continued development of infrastructure in the tundra. The challenge of complex, intertwined natural, social, and built environments on Yamal exemplifies why convergent science is necessary to tackle associated research questions.

3.4 Social System

Within this rapidly changing natural environment on Yamal is a complex social system including indigenous Nenets families living as nomadic reindeer herders and semi-nomadic fishermen, small majority-indigenous villages, and shift workers at infrastructure facilities. The number of permanent residents in Yamal is ca. 17,000 people, and almost 13,000 of them are indigenous peoples (mainly Nenets). The official number of shift workers is 25,000, but we assume that this data is underestimated (Loginov et al., 2020). Approximately 5,500 indigenous people engage in reindeer herding and fishing in the tundra, in a fully nomadic lifestyle with yearly migrations on reindeer sledges of up to 1,200 km (Terekhina & Volkovitsky, 2023). Yamal is one of the few places on the planet where people maintain the kind of nomadic pastoralism where moving is the norm rather than the exception. From the 1960s to present, the domestic reindeer population has grown from between 103,100 and 175,300 in 1990 (Makeev et al., 2014) to approximately 225,000 today (Terekhina & Volkovitsky, 2023). Previous studies (Forbes et al., 2009; Stammler, 2005) argued for strong resilience of reindeer herding lifeways within Yamal socio-ecological systems.

Today, independent households privately manage 80–90% of domestic reindeer in Yamal, while the rest of the herds belongs to a municipal enterprise (one former collective soviet “sovkhoz” that still remained in 2021). While the herding families spend most of their time migrating, most tundra people are registered in one of the villages of Yamalskiy district and children attend school there: Yar-Sale (the district center), Salemal, Syunai-Sale, Panayevsk, Novyi Port, Mys Kamenniy, and Seyakha (Figure S1c in Supporting Information S1). These villages contain core social institutions and infrastructure including administrators, schools, and clinics (Stammler, 2005; Ziker, 2002). Villagers maintain social and cultural relations with their nomadic relatives (Liarskaya, 2016; Volzhanina, 2017). Larger municipalities, such as Salekhard (the capital of YaNAO), have more complex social organization and infrastructure and indigenous leaders (natsional'naia intelligentsia) often reside there.

Beyond the growing development of infrastructure in the tundra, industrialization of Yamal has affected indigenous peoples living in Yamal in other profound ways such as through improved connectivity (expanded cellular network coverage (Stammler, 2009; Stammler, 2016), abundance of government sponsored satellite phones provided for the tundra families, transportation options and government sponsored train tickets) as well as expanded options for goods delivery, fuel supply, healthcare, veterinary services, and selling reindeer meat and fish to shift workers, etc.

Rapid climate change is increasingly impacting reindeer herders (Stammler & Ivanova, 2020). One type of event, rain-on-snow, leads to icing on pastures and inaccessibility of forage for reindeer in winter (Bartsch et al., 2010; Forbes et al., 2016; Sokolov, Sokolova, et al., 2016; Volkovitskiy & Terekhina, 2020). The most prominent of these entailed a loss of approximately 60,000 reindeer on the south Yamal Peninsula during the 2013–2014 winter. In 2018/2019, herds experienced local icings. In the winter of 2020–2021, up to 15,000 domestic and an unknown number of wild reindeer perished in northern Yamal, despite supplemental feed delivered to tundra by the regional authorities using specially chartered aircrafts. A convergent science approach lends itself to understanding the trade offs for various strategies that Nenets reindeer herders employ for dealing with these increasing climatic challenges.

2.1 Satellite Remote Sensing Data

We used six long-term and three short-term satellite remote sensing datasets in this study, covering four types of vegetation indicators - NDVI, LAI, SIF and VOD (Table 1). These remote sensing products include: (a) VIP15 NDVI product (1981–2014), which has 0.05° and 15-day resolutions and is developed by harmonizing the observations of Advanced Very High Resolution Radiometer (AVHRR) from 1981 to 1999 and Moderate Resolution Imaging Spectroradiometer (MODIS) C5 from 2000 to 2014 (Didan et al., 2015); (b) GIMMS NDVI3g (1981–2015), which has 15-day and 1/12° resolutions, and was produced by aggregating daily AVHRR surface reflectance (Pinzon & Tucker, 2014); (c) GIMMS LAI3g product (1981–2015), which was further produced by GIMMS NDVI3g product using a neural network algorithm (Zhu et al., 2013); (d) PKU GIMMS NDVI (1982–2020), which is a new version of GIMMS NDVI product produced by a machine learning model incorporating Landsat images (Li et al., 2023); (e) GLASS LAI product (1981–2018), which has 8-day temporal resolution and 0.05° spatial resolution, and was reconstructed by combing AVHRR LAI from 1981 to 1999 and MODIS LAI from 2000 to 2018 using a bidirectional long short-term memory (Bi-LSTM) model (Ma & Liang, 2022); (f) GLOBMAP LAI (1982–2019) dataset at a spatial resolution of ∼0.07°, covering the period from 1982 to 2019, has half-month (1982–2000) and 8-day (2001–2019) temporal resolutions, and was produced by establishing a pixel-level AVHRR Simple Ratio (SR)-MODIS LAI relationship (Liu et al., 2012); (g) MOD13C1 NDVI (2000–2020) product (C61), which has a 0.05° spatial resolution and a 16-day temporal resolution (Didan & Munoz, 2019); (h) OCO2 SIF product (2000–2020) at resolutions of 0.05° and 4 days, which was generated from MODIS surface reflectance and OCO2 SIF data using a neural network approach (Zhang et al., 2018); and (i) VOD (1987–2017) dataset, which was produced by merging observations from multiple microwave sensors at daily temporal resolution and 0.25° spatial resolution (Moesinger et al., 2020). We standardized these satellite remote sensing products to a 0.5° spatial resolution by using pixel aggregation (PA) method and to a monthly temporal interval by using maximum value composite (MVC) method (Ma et al., 2022; Tian et al., 2021). Considering the high similarity among GIMMS-version datasets (Figure S1 in Supporting Information S1), we only used GIMMS NDVI3g, along with other three independent long-term satellite remote sensing products (i.e., VIP15 NDVI, GLASS LAI and GLOBMAP LAI) for the main analysis.

2.2 LAI Estimates

We used monthly LAI estimates of 12 DGVMs from TRENDY v9 (Friedlingstein et al., 2022; Sitch et al., 2015) to check the consistency of the changes in the vegetation IAV derived from satellite remote sensing products and model measurements. The 12 DGVMs include CABLE, CLASSIC, CLM5, ISAM, ISBA, JULES, LPJ, LPX, ORCHIDEE, SDGVM, VISIT and YIBs (Table S1 in Supporting Information S1), which were driven by monthly CRU or 6-hourly CRU-JRA gridded climate datasets, as well as dynamic atmospheric CO₂ concentrations. These models provide LAI estimates under four different scenarios, namely no change (S0), varying CO₂ only (S1), varying CO₂ and climate (S2), and varying CO₂, climate and land cover change (S3). We analyzed CV trends from LAI estimates of S3 scenario in this study. To match the period of data availability with satellite remote sensing observations, we used LAI estimates from 1982 to 2014.

2.3 Quantifying the Vegetation IAV Changes

We used the coefficient of variation (i.e., CV), the ratio of the standard deviation to the mean, to indicate vegetation IAV for each dataset. The use of CV is meant to remove the differences in magnitude between datasets (i.e., NDVI, LAI, SIF and VOD) to ensure the comparability among the results. To further quantify the CV changes, we first defined the growing season as the period when the mean daily air temperature is above zero (Jiang et al., 2017; Smith et al., 2019) for each vegetated pixel (Figure S2 in Supporting Information S1). We adopted a “methodological growing season” (Körner et al., 2023) as it allows for easy extraction and reproducibility without incurring uncertainties from more complex definitions of phenology. Meanwhile, this definition effectively eliminates the frozen period when there is no vegetation growth, and aligns well with the phenological dates extracted by MODIS NDVI (Leeper et al., 2021). We obtained yearly composites by summing all values within the growing season (Piao et al., 2020; Zhu et al., 2016). Subsequently, we calculated the CV for every 10-year moving window and assigned the CV value to the middle year of the window. The Theil-Sen method was used to estimate the CV trend (i.e., slope) (Wang et al., 2019), and its significance (i.e., p value) was determined using a two-tailed Student’s t-test (Jiang et al., 2019; Xu et al., 2022). Lastly, the non-parametric Mann-Kendall test was used to detect whether a significant monotonic increasing or decreasing trend exists (Jiang et al., 2019).

2.4 Classifying the Levels of Consistency of Remote Sensing-Based CV Trends

Considering potential inconsistency of CV trends across four long-term satellite remote sensing products (i.e., VIP15 NDVI, GIMMS NDVI3g, GLASS LAI and GLOBMAP LAI), we classified the levels of consistency of CV trends for each vegetated pixel using the following criterion (Kause et al., 2022; Xu et al., 2022): (a) “virtually certain” (CER) means the signs of CV trends derived from all four long-term satellite products are the same and significant (Mann-Kendall test, p ≤ 0.1); (b) “likely” (LIK) if the signs of the CV trends are the same and significant in any three satellite products; (c) “about likely as not” (ALN) if the signs of the CV trends are the same and significant in any two satellite products; (d) “possibly” (POS) if only one satellite product yields significant CV trend, and others are insignificant; (e) “no change” (NOC) if no significant CV changes were detected in all four satellite products; and (f) “conflicting” (CON) if the observed CV trends are conflicting with each other (i.e., significant positive and negative CV trends were detected simultaneously across the four products). We further assigned “+” and “−” signs to the consistency levels that refer to the direction of CV trends (i.e., positive and negative trends, respectively). Overall, the levels of consistency were classified into four types, positive, that is, CER (+), LIK (+), ALN (+) and POS (+), negative, that is, CER (−), LIK (−), ALN (−) and POS (−), no change (NOC), and conflicting (CON), respectively. Based on previous studies, NDVI and LAI are both indicators of canopy structure and greenness, and they are strongly related (Wang et al., 2022; Zeng, Hao, et al., 2022). Therefore, we deem that there is a clear physiological and statistical basis to treat them equally. In addition, we used both indicators to augment the number of available long-term datasets for our analysis, as some datasets (e.g., VIP15) only contain NDVI, while others (e.g., GLOBMAP) only include LAI.

2.5 Identifying Factors Driving the IAV Changes

Based on the findings from previous studies (Baldocchi et al., 2016; Berdugo et al., 2022; Forkel et al., 2016; Gray et al., 2014; Luo & Keenan, 2022; Piao et al., 2015; Zhu et al., 2016), we investigated several factors that may serve as potential drivers for the IAV changes. These factors are aridity index (AI), mean annual air temperature (MAT), mean annual precipitation (MAP), land use and land cover change (LUCC), mean monthly LAI, the trend of LAI (LAI_trend) and nitrogen deposition (N deposition). AI was obtained from the Global Aridity Index dataset at a 30 arc-second resolution from 1970 to 2000 (Zomer et al., 2022), in which lower AI values mean drier conditions. MAT and MAP data were extracted from the gridded CRU-JRA V2.1 dataset (6-hr and 0.5° resolutions) (Harris et al., 2020). LUCC intensity was obtained using a 300-m ESA CCI land cover type product with 22 land cover classes from 1992 to 2019. We first determined whether land cover change occurred in each 300-m pixel from 1992 to 2019, and then quantified the LUCC intensity as the proportions of land cover change in each 0.5° spatial grid. The mean annual LAI was the average of GLASS LAI and GLOBMAP LAI from 1982 to 2014, and the LAI_trend (i.e., slope) was estimated using Mann-Kendall and Theil-Sen method for each vegetated pixel. Most global vegetated regions show significant greening trends (i.e., LAI_trend > 0) over the past decades (Chen, Chi, et al., 2019; Zhu et al., 2016). This implies that terrestrial ecosystems have accumulated more green leaves for carbon fluxes and in response to climate change, which may induce IAV changes. The N deposition was the average of the nitrogen deposition dataset (2° × 2.5° grid resolution) from 1984 to 2016 (Ackerman et al., 2019). All the above datasets were resampled into a 0.5° spatial grid.

We first assessed the significance of each factor in driving the vegetation IAV by using a two-sample t-test method, in which we examined whether the high and low values of each factor caused statistically different CV trends. For each factor, we used different criteria for the binary classification of high and low values, that is, the 50% for LUCC, the zero for LAI_trend, and the global average for other drivers. After that, we selected the four factors (i.e., AI, MAT, LAI_trend and N deposition; see Results) that have significant impacts on the CV trend and quantified their respective impacts using a multiple regression method: