Impacts of abiotic and biotic factors on tundra productivity near Utqiaġvik, Alaska

Earlier snowmelt, warmer temperatures and herbivory are among the factors that influence high-latitude tundra productivity near the town of Utqiaġvik in northern Alaska. However, our understanding of the potential interactions between these factors is limited. MODIS observations provide cover fractions of vegetation, snow, standing water, and soil, and fractional absorption of photosynthetically active radiation by canopy chlorophyll (fAPARchl) per pixel. Here, we evaluated a recent time-period (2001–2014) that the tundra experienced large interannual variability in vegetation productivity metrics (i.e. fAPARchl and APARchl), which was explainable by both abiotic and biotic factors. We found earlier snowmelt to increase soil and vegetation cover, and productivity in June, while warmer temperatures significantly increased monthly productivity. However, abiotic factors failed to explain stark decreases in productivity during August of 2008, which coincided with a severe lemming outbreak. MODIS observations found this tundra ecosystem to completely recover two years later, resulting in elevated productivity. This study highlights the potential roles of both climate and herbivory in modulating the interannual variability of remotely retrieved plant productivity metrics in Arctic coastal tundra ecosystems.


Introduction
Monitoring, quantifying, and understanding changes in vegetation, snow cover and permafrost thaw in high latitude regions is fundamental for predicting how these regions and their inhabitants are impacted by climate change and how humans will need to adapt (Hinzman et al 2005, Carroll et al 2011, Arp et al 2013, Vincent et al 2013).Warming of northern latitudes has coincided with shifts in tundra vegetation (Mekonnen et al 2021), disappearing ponds and lakes (Andresen andLougheed 2015, Jones et al 2020), earlier snowmelt, lengthening of the growing season and start of spring greening and increase in tundra productivity (Myneni et al 1997, Wu et al 2022, Liu et al 2023).
High-latitude ecosystems are characterized by a short growing season and a long snow season (Serreze 2010, Zhang 2021).Warming and changes in the physio-hydrological environment has not only altered the timing of snowmelt and increased permafrost active layer thickness, but also plant production.The normalized difference vegetation index (NDVI) products from the advanced very high-resolution radiometer (AVHRR) have been often used to assess high latitude greening/browning trends.These studies used NDVI as a proxy of tundra productivity, and found pan-Arctic trends in NDVI to increase with the summer warmth index (SWI; the sum of average monthly temperature greater than 0 Though knowledge of the spatiotemporal linkages between abiotic factors and tundra productivity is growing, the short and long-term impacts of acute stocastic biotic disturbances is limited.In northern Alaska near the town of Utqiaġvik, the impact of the dominant herbivore (i.e.brown lemming; Lemmus trimucronatus) on vegetation productivity was intensively studied during the early 1950s by field ecologists (Norton 2001).These early plot level experimental field studies reported that the digging of rhizomes by lemmings over the winter dramatically reduced the standing crop of forage by 20%-50% during the following summer (Thompson 1955, Schultz 1964).Field surveys during the recent 2008 lemming population boom (Villarreal et al 2012) supported these early observations, dramatically reducing vegetation productivity metrics, GPP and NDVI (Lara et al 2017, Pitelka andBatzli 2007).Though short-term responses of defoliation directly reduced productivity, the long-term impacts of such intense biotic disturbances are more elusive and difficult to predict (Lara et al 2017).MODIS and AVHRR pixels in this region are often mixtures of soil, ponds, vegetation, and/or snow/melt.Soil, ponds, and snow/melt have confounding impacts on MODIS and AVHRR NDVI.In addition, NDVI has limited capability to distinguish between photosynthetic vegetation and nonphotosynthetic vegetation (Zhang et al 2020(Zhang et al , 2021)).
Understanding the complex interactions between abiotic and biotic factors within data-rich subregions of the Arctic may improve our ability to interpret decadal-scale patterns and trends in remotely retrieved productivity metrics spanning tundra landscapes.Here, we evaluate the seasonal to decadal-scale variability in MODIS derived productivity metrics between 2001 to 2014, during a period of extraordinarily high variability in climate and environmental conditions (i.e.air temperature, downwelling and upwelling irradiance, and snowmelt date) and within the last known lemming population outbreak near Utqiaġvik (i.e. 2008).The productivity metrics used in this study include the fractional absorption of photosynthetically active radiation (PAR) by canopy chlorophyll (fAPAR chl ) and the absorption of PAR (APAR) by canopy chlorophyll (APAR chl ) (Zhang 2021; table 1).The fAPAR chl is superior to NDVI in GPP estimation (Zhang et al 2005).Specifically, we address the following questions: (1) how did earlier snowmelt and air temperature affect the Utqiaġvik (formerly Barrow) tundra productivity and surface cover fractions?and (2) how did the 2008 lemming outbreak impact remotely retrieved vegetation productivity (e.g.fAPAR chl and APAR chl )?This study advances our ability to interpret spatial and temporal patterns and trends in remotely retrieved tundra productivity metrics.

Meteorological measurements
Our study examined the terrestrial land surface area of a 7 × 7 km 2 region surrounding the US-Brw site (71.32 • N, 156.61 • W) near the town of Utqiaġvik on the northern Barrow Peninsula of Alaska.This region is dominated by thermokarst lakes and drained thaw lake basins, which are composed of polygonal ice-wedge tundra (Brown et al 1980, Lara et al 2015, 2018b) (figure 1).NOAA's Earth System Research Laboratory (ESRL) at Utqiaġvik measures downwelling and upwelling irradiance, wind, precipitation, and air temperature (Longenecker 2017, Vasel 2019).From 2001 to 2014, this area had a mean precipitation of 44.5 mm as rain during a three-month period from July through September.Mean annual temperature was −10.4 • C, with lowest temperature in February (−24.6 • C) and highest temperature in August (+4.2• C).The snow season spans from October to mid-May of the next year.Melting of snow increases soil moisture due to continuous permafrost and the low topographic relief (Hobbie 1984, Jepsen et al 2013, Carroll and Loboda 2017, Nitze et al 2017).
Upwelling irradiance increases with downwelling irradiance due to the presence of snow during the winter and early spring, then declines quickly with snowmelt.We utilized the method developed by Stone et al (2002) and the NOAA ESRL measurements of downwelling and upwelling irradiance to determine the snowmelt date: the first day of year (DOY) when the albedo started to be less than 0.3.The NOAA ESRL measurements of air temperature were used to compute SWI (Jia et al 2003;table 1).
We leveraged site-level observations near Utqiaġvik to identify the timing of high lemming population densities (Johnson et al 2011, Villarreal et al 2012, Lara et al 2017).Lara et al (2017) reviewed the literature to chronologically record all known lemming outbreaks at our site, which occurred in 1946, 1949, 1953, 1956, 1960, 1965, 1971, 1976, 1981, 1986, 1991, 1996, 2000, 2003, and 2008.Due to the logistical challenges of continuous trapping needed to estimate annual lemmings ha −1 in remote tundra regions, the long history of outbreak years remains among the most reliable and longest running boombust observations in the Arctic.We leverage these ground-based observations to advance our ability to remotely detect lemming outbreaks from space.

Satellite observations
At high latitudes, there are frequent MODIS overpasses that increase opportunities for clear sky data collections.MODIS data of land bands 1-7 at a 500 m resolution (Masuoka 2019) from 2001 to 2014 were employed for the study region: blue (459-479 nm), green (545-565 nm), red (620-670 nm), near infrared (NIR 1 : 841-875 nm; NIR 2 : 1230-1250 nm), and shortwave infrared (SWIR 1 : 1628-1652 nm; SWIR 2 : 2105-2155 nm) (https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/6/MCD19A1).We used high quality surface reflectance processed with the multi-angle implementation of atmospheric correction (MAIAC) algorithm (Lyapustin et al 2018).MAIAC is an advanced algorithm that uses time series analysis and a combination of pixel-based and imagebased processing to improve accuracy of cloud/snow detection, aerosol retrievals, and atmospheric correction by incorporating the bidirectional reflectance distribution function model of the surface.The 500 m surface reflectance data were used to monitor the spatial and temporal dynamics of fAPAR chl of the tundra ecosystem.The fAPAR chl and fractional APAR by canopy non-chlorophyll components (fAPAR non-chl ) were retrieved with surface reflectance of the seven bands and the coupled leaf-vegetationsoil-snow-standing water model (LVS3) (Zhang et al 2020(Zhang et al , 2021)).APAR chl is the product of fAPAR chl and PAR (APAR chl = fAPAR chl × PAR).In addition to the retrieval of the temporal dynamics of these productivity metrics, we also evaluated the coincident patterns in the cover fractions of vegetation (VGCF), soil (SOILCF), standing water (WaterBodyCF), and snow (SNOWCF) derived from MODIS (table 1).The Utqiaġvik landscapes are often characterized by some portion of lakes, ponds, and seasonal freestanding water.The LVS3 model is a comprehensive model that depicts the interactions between leaf, vegetation, soil, snow, and surface water.It originated from the PROSAIL2 model (Zhang et al 2005, 2006, 2009, 2014a, 2014b, 2015, Cheng et al 2014).PROSAIL2 (PROSPECT + SAIL2) is a radiative transfer model that simulates the reflection and transmission of solar radiation by vegetation and soil.The LVS3 model is a useful tool for studying the interactions between different land surface components, and for estimating various vegetation, water, snow and soil metrics using remote sensing data.Monthly VGCF, SOILCF, WaterBodyCF, SNOWCF, fAPAR chl and fAPAR non-chl products in April, May, June, July, August and September were produced from the MODIS observations with the LVS3 model.The chlorophyll molecules in vegetation absorb sunlight for energy to drive photosynthesis, which is a temperature-sensitive process (Raich et al 1991, Zhang et al 2005).Downwelling irradiance reached its peak in June.APAR chl is computed with MODIS fAPAR chl and the NOAA ESRL measurements-based PAR.50% of shortwave radiation is in the PAR spectral interval (Pinker and Laszlo 1992, Frouin and Pinker 1995, Pinker et al 2010;www2.atmos.umd.edu/∼srb/par/webpar.htm).
In additon to MODIS data, we acquired the Landsat TM images for this region during 2001-2014 at a 30 m resolution from Google Earth Engine (Hu et al 2021).July and August fAPAR chl were retrieved.A summertime growth index was further derived as the fAPAR chl difference between July and August.The deviation of the summertime growth in the year 2008 relative to the 14 years' average of the growth index was then computed as a measure to quantify fine-scale patterns of herbivory disturbances.

Cover fractions of vegetation, soil, water and snow
Evaluating interannual patterns of biophysical variables (e.g.vegetation, soil, surface water and snow cover) provides insights into short and long-term patterns of ecological change (figure 3).The 2001-2014 average monthly PAR value declined from June to September.The 14 year average July and August temperature values were higher than the June and September values.Cover fractions of vegetation, soil, surface water and snow generally fluctuated from year to year, but we identified positive and negative trends in vegetation and soil cover fractions in July and August, respectively (figures 3(C) and (D)).

fAPAR chl
Plant chlorophyll content is an important indicator of vegetation productivity in terrestrial ecosystems (Zhang et al 2005(Zhang et al , 2014a)).Typical Arctic tundra annual phenological patterns of fAPAR chl begin with low values following snowmelt, steadily increase from June to August as vegetation growth expands spatially, and decline during the senescence period (figure 4(A)).The long-term average monthly fAPAR chl values in June, July, August and September were 0.03, 0.16, 0.21 and 0.11, respectively.The annual percent changes of June, July, August and September fAPAR chl were 0.55, 1.66, 0.98, and −1.70, respectively (figure 4(B)).

APAR chl
The importance of interannual changes in fAPAR chl (figure 4(A)) to productivity was modulated by the strong seasonal variability in PAR (figure 2(A)).The 14 year averages for both fAPAR chl and PAR in June and September differed greatly, with the fAPAR chl value in September was about 4.2 times the fAPAR chl value in June while the PAR value in June was about 4.5 times the PAR value in September.Therefore, the long-term averages for APAR chl in June and September were similar (figure 4(A)).
The 14 year average monthly APAR chl values in June, July, August and September were 110.30W m −2 , 567.82 W m −2 , 429.96 W m −2 , and 113.10 W m −2 , respectively.Thus, percentages of monthly APAR chl values to the entire growing season (June-September) APAR chl value were 9%, 46.5%, 35.2% and 9.3%, respectively.Interannual variability in APAR chl illustrates the effects of both chlorophyll content and cloudiness on productivity (figure 4(A)).The seasonal curve for APAR chl differed from that of   fAPAR chl , with the APAR chl having sharper, higher, and earlier peak (figure 4(A)).
Since APAR chl is closely related to productivity (Zhang et al 2014a), we show the multiyear variability in APAR chl for June, July, August, September, July-August and the growing season (June-September) (figure 4(C)).There was large interannual variability in the growing season APAR chl , which was mainly driven by variability in the APAR chl for July-August.The large interannual variability previously shown for fAPAR chl in September (figure 4(B)) was dampened for APAR chl due to the low levels of incident PAR in September (figure 4(C)).The annual percent changes of June, July, August and September APAR chl were 1.18, 1.46, 1.34, and −0.29, respectively.The annual percent change of July-August APAR chl was 1.41.The annual percent change of June-September APAR chl was 1.23 (figure 4(C)).

Influence of early snowmelt and temperature on vegetation
Relationships between the snowmelt date and (1) monthly fAPAR chl in June, July, August, and September; (2) monthly APAR chl in June, July, August and September; (3) total July-August APAR chl ; and (4) total June-September APAR chl were evaluated (figures 5(A) and (B)).The correlation between snowmelt date and monthly fAPAR chl in June was significant (figure 5(A): F = 10.01,p = 0.008 < 0.05).The correlation between snowmelt date and monthly APAR chl in June was significant (figure 5(B): F = 6.03, p = 0.030 < 0.05).
Temperature was an essential control of productivity on the Barrow Peninsula (Lara et al 2018a).SWI significantly impacted monthly fAPAR chl for August (figure 6(A): F = 4.84, p = 0.048 < 0.05) and monthly APAR chl for August (F = 13.03,p = 0.004 < 0.05) (figure 6(B)).SWI also had a significant impact on APAR chl for July-August (F = 6.00, p = 0.031 < 0.05).As the contributions of APAR chl in June and September to the growing season total APAR chl were low, SWI had a significant impact on the growing season APAR chl (June-September) (F = 5.30, p = 0.040 < 0.05).

Influence of lemming herbivory on vegetation
Typically, fAPAR chl increased from July to August as observed over the 14 year observation record (figure 4 Value of fAPAR chl in July of 2008 was even greater than that in August due to the severe impact of the 2008 lemming population outbreak on foliage thus chlorophyll contents and fAPAR chl .Figure 8 shows the

Discussion
Over a decade of MODIS observations near Utqiaġvik finds that this tundra region experienced large interannual variability across all vegetation productivity metrics, driven by both abiotic and biotic factors.Due to the large interannual variability in all vegetation metrics in June, July and August, their trends were typically not within the 95% confidence level.Monthly fAPAR chl and APAR chl in June, July and August increased with temperature (figures 6(C) and (D)) and also increased over time (figures 4(B) and (C)), suggesting that warming may increase tundra greening and vegetation productivity.Furthermore, similar to Zona et al (2022), our results of fAPAR chl and APAR chl suggest that an earlier snowmelt may indicate an earlier green-up, and overall greater vegetation productivity, specifically in June and July (figures 5(A) and (B)).These findings support the hypothesis that earlier snowmelt is indeed a predictor of higher Arctic production in June and July (Goetz et al 2005, Rosa et al 2015).
The NDVI from MODIS were computed and described in the supplementary material (figure S1).We found that an earlier snowmelt had a significant impact on monthly NDVI for June (F = 12.74, p = 0.004 < 0.05) (see figure S2(A)), not only due to the impact of earlier snowmelt on vegetation cover fraction (VGCF) for June (F = 9.44, p = 0.010 < 0.05) but also due to its impact on soil cover fraction (SOILCF) (F = 45.86,p = 0.000 < 0.05) (figure 5(C)).Seasonal VGCF, fAPAR chl and APAR chl profiles are different from seasonal NDVI profile (figures 3(A), 4(A) and S1(A)), hinting time series of VGCF, fAPAR chl and APAR chl may potentially provide alternative information for satellite-derived land surface phenology studies (Ye et al 2022, Zhang et al 2022), and for disturbance studies (Zhao et al 2018).
Although interannual controls on vegetation productivity are influenced by the timing of snowmelt, SWI and monthly temperature, our analysis finds local-scale dynamic disturbances confound decadal productivity trends.Due to the impact that lemmings can have on vegetation composition, nutrient cycling, and food web dynamics (Batzli et al 1980, Ims andFuglei 2005)  Following the lemming outbreak, MODIS observations show a gradual recovery in 2009, and a complete recovery by 2010, resulting in a 'greener' landscape with greater fAPAR chl than 2007 (figure 7(A)).This two-year vegetation recovery period (figure 7) is supported by field observations (Johnson et al 2011, Villarreal et al 2012, Lara et al 2017) and predicted by the nutrient recovery hypotheses (Pitelka 1964, Batzli et al 1980, Pitelka and Batzli 2007).This hypothesis describes the process of rodent boom-bust dynamics, whereby high rodent population densities (i.e.boom) decimate the vegetation, increasing bare ground and altering surface energy balance (increasing thaw depth), which alters nutrient dynamics, resulting in reduced quality and quantity of forage.Lemming populations crash (i.e.bust) and gradually recover only when plant production increases, restoring surface energy balance and nutrient dynamics to levels that can meet the nutritional demands of high lemming population densities (Pitelka and Batzli 2007).However, perhaps most importantly, results find the lemming outbreak year to be the only year on record where monthly fAPAR chl did not increase from July to August (figures 7(A) and (B)).The rare stability of monthly productivity metrics presents a new and exciting opportunity to use remote sensing to detect lemming outbreaks across vast regions of the Arctic.We will explore in this direction in the future.

Conclusions
We determined the importance of both abiotic and biotic factors in interpreting over a decade of trends, interannual variability and seasonal phenological patterns in vegetation productivity across the Barrow Peninsula.We found productivity metrics to increase with temperature and the SWI on decadal, interannual, and seasonal time scales, though earlier snowmelt consistently resulted in earlier green-up and overall greater annual vegetation productivity.However, these traditionally important drivers/indicators can be shaken up during disturbance years, as highlighted during the 2008 lemming outbreak.In line with field observations, experiments, and the nutrient recovery hypothesis, we were able to characterize the vegetation recovery period with fAPAR chl following disturbance.This satellite remote sensing analysis with LVS3 provided useful information for documenting surface changes in response to top-down and bottom-up driving factors, which may provide scientific support for understanding and managing high latitude ecosystems.

Figure 2 .
Figure 2.For the years 2001-2014: (A) the monthly PAR (W m −2 d −1 ) for June, July, August and September; (B) the snowmelt dates in day of year (DOY) determined with the method developed by Stone et al (2002) and the NOAA ESRL measurements of upwelling and downwelling irradiance; (C) the monthly temperature ( o C) for June, July, August and September; and (D) the summer warmth index (SWI, unit: o C).

Figure 3 .
Figure 3. Cover fraction dynamics of vegetation, soil, water and snow: (A) average seasonal VGCF, SOILCF, WaterBodyCF and SNOWCF in 2001-2014 for the study area; (B) time series of VGCF, SOILCF and WaterBodyCF in June of 2001-2014; (C) time series of VGCF, SOILCF and WaterBodyCF in July; (D) time series of VGCF, SOILCF and WaterBodyCF in August; and (E) time series of VGCF, SOILCF and WaterBodyCF in September.

Figure 5 .
Figure 5. Relationships between snowmelt date in day of year (DOY) and (A) fAPAR chl in June, July, August and September; (B) APAR chl in June, July, August and September; and (C) VGCF, SOILCF, SNOWCF and WaterBodyCF in June.
(A)).During the summer of 2008, a relatively high lemming population outbreak caused severe vegetation defoliation (Johnson et al 2011, Villarreal et al 2012, Lara et al 2017).fAPAR chl failed to increase from July to August in 2008 (figures 7(A) and (B)).

Figure 6 .
Figure6.Impacts of SWI and monthly temperature on fAPAR chl and APAR chl : (A) relationships between SWI (unit: o C) and monthly fAPAR chl in June, July, August, September; (B) relationships between SWI and total APAR chl for June, July, August, September; and (C) relationships between monthly temperature (unit: o C) and monthly fAPAR chl in June, July, August, September; and (D) relationships between monthly temperature and total APAR chl for June, July, August, September.
, the severe lemming outbreak of 2008 was indeed noticeable from MODIS observations (figure 7).During the winter of 2007-2008, lemming likely grew rapidly under the snow and their population exploded during the summer of 2008 (Johnson et al 2011, Lara et al 2017).A high proportion of

Figure 8 .
Figure 8. Spatial pattern of the anomaly in summertime vegetation growth for the year 2008 relative to the long-term average over 2001-2014.Summertime vegetation growth was quantified as the difference in fAPAR chl between August and July (i.e.fAPAR chl _Aug-fAPAR chl _Jul).Warm colors indicate negative anomalies, attributable largely to the lemming outbreak.

Table 1 .
Acronyms used in this study.