Attribution of the spatial heterogeneity of Arctic surface albedo feedback to the dynamics of vegetation, snow and soil properties and their interactions

The Arctic warming rate is triple the global average, which is partially caused by surface albedo feedback (SAF). Understanding the varying pattern of SAF and the mechanisms is therefore critical for predicting future Arctic climate under anthropogenic warming. To date, however, how the spatial pattern of seasonal SAF is influenced by various land surface factors remains unclear. Here, we aim to quantify the strengths of seasonal SAF across the Arctic and to attribute its spatial heterogeneity to the dynamics of vegetation, snow and soil as well as their interactions. The results show a large positive SAF above −5% K−1 across Baffin Island in January and eastern Yakutia in June, while a large negative SAF beyond 5% K−1 is observed in Canada, Chukotka and low latitudes of Greenland in January and Nunavut, Baffin Island and Krasnoyarsk Krai in July. Overall, a great spatial heterogeneity of Arctic land warming induced by positive SAF is found with a coefficient of variation (CV) larger than 61.5%, and the largest spatial difference is detected in wintertime with a CV > 643.9%. Based on the optimal parameter-based geographic detector model, the impacts of snow cover fraction (SCF), land cover type (LC), normalized difference vegetation index (NDVI), soil water content (SW), soil substrate chemistry (SC) and soil type (ST) on the spatial pattern of positive SAF are quantified. The rank of determinant power is SCF > LC > NDVI > SW > SC > ST, which indicates that the spatial patterns of snow cover, land cover and vegetation coverage dominate the spatial heterogeneity of positive SAF in the Arctic. The interactions between SCF, LC and SW exert further influences on the spatial pattern of positive SAF in March, June and July. This work could provide a deeper understanding of how various land factors contribute to the spatial heterogeneity of Arctic land warming at the annual cycle.


Introduction
Surface albedo refers to the ratio of the total reflected radiation to the incident radiation for a certain surface (Henderson-Sellers and Wilson 1983), and how its variations can significantly regulate regional temperature anomalies (Dai et al 2019, Alessandri et al 2021. Generally, surface albedo feedback (SAF) includes positive and negative effects. A negative SAF refers to an increase in highly reflective factors (e.g. snow and sea ice) of the underlying surface and a decrease in the absorption of incoming solar radiation, thus strengthening surface cooling (Hall andQu 2006, Alessandri et al 2021). In contrast, the phenomenon of increased temperature triggered by reduced albedo is often referred to as a positive SAF (Hall and Qu 2006). It has long been known that the loss of snow and ice under global warming greatly contributes to the decrease in Arctic surface albedo (Walsh 2014, Andry et al 2017, Huang et al 2017, Jakobs et al 2020, which enhances the absorption of solar radiation and thus amplifies the warming rate over the Arctic (see figure S1 (available online at stacks.iop.org/ERL/17/014036/mmedia), feedback loop). Therefore, understanding the mechanisms behind historical positive SAF is fundamental for predicting future climate warming in the Arctic due to global climate change.
Previous studies have mainly focused on quantifying the strengths of SAF and its components, as well as its seasonal cycle at various time and space scales (Hall and Qu 2006, Wegmann et al 2018, Thackeray and Hall 2019, Alessandri et al 2021. For example, Qu and Hall (2014) indicated that approximately 5%-10% of the global warming rate is contributed by the variation in SAF via changes to the top of atmosphere radiation balance. However, SAF can exert a greater impact on regional climate and hydrological cycles because this feedback tends to be larger in the regions where snow and ice dynamics are susceptible to climate change (Bowman et al 2018, Kim et al 2018, Thackeray et al 2018, Webb et al 2021. A number of studies have estimated that an approximately 1% reduction in land albedo over the Northern Hemisphere could cause an average of 1 • C of near surface warming, equivalent to a radiative effect of approximately 0.1 W m −2 K −1 (Qu and Hall 2006, 2014, Flanner et al 2011, Fletcher et al 2014. As would be expected due to low snow and ice cover, Qu and Hall (2014) found a low SAF (0.1% K −1 -0.3% K −1 ) during summertime compared to other seasons (0.6% K −1 -0.9% K −1 ) across Northern Hemisphere extratropical land areas. Similarly, over the Arctic, the SAF exerts an impact mainly in spring and/or summer, which is closely associated with the seasonal cycling of sea ice area (Andry et al 2017). Spatially, both observational and modeling analyses reveal stronger feedback strengths in high latitudes from 50 • N to 70 • N (Fletcher et al 2012, 2014, Thackeray et al 2021 with a large positive SAF (+0.87 W m −2 K −1 ) in snow-dominated areas (annual mean areal fraction of snow cover >10%) (Alessandri et al 2021).
In recent decades, a considerable reduction in snow cover extent/area has been detected in spring and summer by satellite-derived observations in polar and high mountainous regions due to climate warming (Estilow et al 2015, Hartfield et al 2018, Notarnicola 2020. Moreover, extensive greening of the Earth was reported under anthropogenic climate warming (Forzieri et al 2017), which could influence SAF in snow-covered regions (e.g. polar areas and high mountains) by altering the masking effect of vegetation over snow and land surface albedo (Loranty et al 2014, Yu et al 2021. In addition to snow and vegetation, albedo can also be modulated by soil parameters, including soil type (ST) and soil moisture. For instance, dry soil often shows high albedo, and vice versa for wet soil (Yang et al 2020). In addition, a smaller soil particle size tends to increase soil albedo because incident radiation is more easily trapped in large soil aggregates with more inter-aggregate spaces and cracks (Cierniewski et al 2013). These results suggest that the dynamics of snow, vegetation and soil properties could interact with each other and modulate the spatial pattern of land surface albedo, affecting the spatial heterogeneity of Arctic land warming induced by positive SAF. To date, however, quantitative attribution of the spatial pattern of Arctic land warming induced by positive SAF across Arctic land areas is still lacking.
To fill this gap, this study aims to quantify the strengths of SAF across the Arctic and to identify the key influencing factors governing its spatial heterogeneity during the historical period of 1982-2015. The three key scientific questions addressed in this study are as follows: (a) how is the seasonal SAF spatially distributed across Arctic land near the surface? (b) How do the dynamics of vegetation, snow and soils control the spatial pattern of seasonal positive SAF? (c) How do vegetation, snow and soils interact with each other and affect the positive SAF? The results can enhance our understanding of the regional differences in SAF across the Arctic and can have large implications for predicting future Arctic climate under global warming.

Study area
Here, the Arctic is defined as the areas north of the Arctic Circle (66 • 34 ′ ), including the Arctic Ocean, adjacent seas, and parts of Alaska, Northwest Territories, Nunavut, Chukotka, eastern Yakutia, western Yakutia, Krasnoyarsk Krai, Yamalo-Nenets Autonomous Okrug (YaNAO) and Greenland (figure 1). Cool summers and cold winters are typical features of the Arctic climate. Generally, most regions of the Arctic receive annual precipitation below 500 mm, mainly in the form of snow (Przybylak 2003

Land surface albedo dataset
The Global LAnd Surface Satellite (GLASS) phase-2 broadband surface albedo product is adopted in this study (www.glass.umd.edu/Download.html), which is characterized by long time series and spatiotemporal continuity (Liang et al 2013, 2021). This product is derived from advance very-high-resolution imaging spectroradiometer (AVHRR) data at a spatial resolution of 0.05 • and a temporal resolution of 8 d. During the process of generating the GLASS phase-2 surface albedo product, the broadband surface albedo (shortwave, visible, and near infrared) is first estimated using the direct-estimation algorithm from the cloud-detected moderate-resolution imaging spectroradiometer and AVHRR data and then filled and fused by statistics-based filter approaches (Liang et al 2013, 2021, Jia et al 2018. The GLASS albedo product provides information on the directional hemispherical reflectance (black sky albedo) and bihemispherical reflectance (white sky albedo). The shortwave broadband (0.3-5.0 µm) black sky albedo is used in this study and is defined as the reflection of direct radiation in the absence of a diffuse component of illumination (Liang et al 2021).

NDVI dataset
The GIMMS3g normalized difference vegetation index (NDVI) is produced by the National Oceanic and Atmospheric Administration of the USA from AVHHRR data (Tucker et al 2005) and is available at https://ecocast.arc.nasa.gov/data/pub/gimms/3g. v1/. The 15 d GIMMS3g NDVI data have a spatial resolution of 1/12 • (∼8 km) covering the period from 1982 to 2015. This set of data has eliminated the effects of sensor replacement and orbital drift during the production process and is widely applied for measuring ecosystem carbon cycles, vegetation phenological changes and climate responses (Kong et al 2016, Pan et al 2018, Yu et al 2021. In this study, the maximum value composite (MVC) technique is used to produce a monthly NDVI. The specific introduction of the MVC method can be seen in the supplementary information (SI).

Land parameter data
The 2 m air temperature (T2m), snow cover fraction (SCF), soil water content (SW) and ST data are obtained from the land component of the 5th generation of European ReAnalysis (ERA5-Land) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/ reanalysis-era5-land?tab=form). ERA5-Land is mainly used to describe the evolution of land parameters in an accordant mode in recent decades (Muñoz-Sabater et al 2021). Furthermore, ERA5-Land shares most parameterizations with ERA5, thus, ensuring the use of a state-of-the-art land surface model (Muñoz-Sabater et al 2021). Another great merit of ERA5-Land is its improved spatial resolution to 9 km from the 31 km ERA5 and 80 km ERA-Interim products (Muñoz-Sabater et al 2021).

Circumpolar Arctic vegetation map (CAVM)
The CAVM is the first vegetation map covering the entire Arctic based on a single and unified legend (www.arcticatlas.org) (CAVM Team 2003). The CAVM is a vector map compiled from handdrawn polygons that are interpreted by a geobotanical method (Walker et al 2005). This map contains 15 vegetation types using a relatively simple legend (Walker et al 2005). In addition, the map also includes the outlines of different bioclimate subzones, lake cover, physiography and substrate chemistry. The land cover types (LCs) (mainly vegetation cover types) and soil substrate chemistry (SC) are used in this study. The spatial patterns of LC and SC are shown in figure S2 in the SI.

Data preprocessing
All the above raster data are reprojected and resampled to the 25 km equal-area scalable Earth grid (EASE-grid), which is a widely used format for data analysis and visualization for polar regions (Atlaskina et al 2015, Yu et al 2021, Yu and Leng 2022. For albedo data, first, the ∼0.05 • resolution datasets are reprojected to the EASE 5 km grid using the nearest neighbor approach such that each 25 km grid contains 5 × 5 smaller grids. Second, the data are aggregated using an averaging method to the 25 km EASE grid in ArcGIS software. Similarly, the T2m, SCF, SW and ST data are homogenized over the Arctic to a 25 km spatial resolution based on the bilinear interpolation method.

Evaluation of the strength of surface albedo feedback (SAF)
Following Hall and Qu (2006) and Qu and Hall (2014), the strength of SAF is defined as the relative percentage change in albedo ( ∆Albedo×100 Albedo ), divided by the corresponding change in near surface (2 m) air temperature (∆T2m). The specific calculation formula is as follows: where ∆ denotes the month-to-month change, which represents 1 month to the next successive month. From this equation, it can be seen that an increase in albedo from 1 month to the next month, as expected in autumn with increasing snow and ice cover, leads to a negative SAF, while a decrease in albedo, as expected in spring, leads to a positive SAF.

Evaluation of the spatial heterogeneity of positive SAF
Following previous studies (e.g. Hall andQu 2006, Qu andHall 2014), negative values of SAF indicate a positive SAF, while a decrease in albedo could enhance the near surface temperature. Consequently, we only evaluate the spatial heterogeneity of positive SAF over the Arctic in this study. The coefficient of variation (CV), which describes the relative deviation of a series of data from its mean value, is adopted to evaluate the spatial heterogeneity of positive SAF in this study. The calculation formula is as follows: where P-SAF SD indicates the standard deviation of positive SAF across the Arctic land areas, and P-SAF M denotes the averages of these positive SAFs.

Changing trend calculation and test
Based on a regression analytical method, the changing trend of CVs of positive SAF can be calculated in different months for the 1982-2015 period. The trends are estimated using equation (3): where i is the serial number of each year; X i is the value of CV in year i; and n is the accumulative number of years during the study periods. The Mann-Kendall trend test, commonly known as the Kendall statistic, is a non-parametric evaluation of monotonic trends of a variable over time that does not require measurements to be normally distributed or the trend to be linear (Mann 1945, Kendall 1948. In the Mann-Kendall test, the null hypothesis H 0 assumes that if no trend exists, the n measurements x 1 , ···, x n obtained over time are realizations of the independently and identically distributed random variables X 1 , ···, X n . The alternative hypothesis H 1 for a two-sided test assumes that the distributions of X k and X j are not independently and identically distributed. For all k, j ⩽ n; with k ̸ = j. The test statistics, with zero mean and asymptotically normal variance, can be calculated using equations (4) and (5): Sign where S is the Gaussian distribution, and the mean value is 0. The variance can be calculated by Var (S) = UW Uni-weakening: impacts of single variables are univariably weakened by the interaction.
BE Bivariate-enhancement: impact of single variables are bivariably enhanced by the interaction.
NE Nonlinear-enhancement: Impacts of variables are nonlinearly enhanced.
. In cases where the sample size n > 10, the standard normal variate Z is computed using: Thus, in a two-sided trend test, H 0 is not rejected if |Z| is smaller than Z α/2 at the α level of significance. A positive Z value indicates an increasing trend, while a negative Z value indicates a decreasing trend.

Optimal parameter-based geographic detector (OPGD) model
The OPGD model is adopted to investigate the spatial heterogeneity of a geographic phenomenon and potential influencing factors (Song et al 2020). This technique is based on the following hypotheses: if the driving factors have a significant influence on the warming of Arctic land, the spatial patterns of the independent and dependent variables should be similar (Song et al 2020). Compared to the traditional linear model, the OPGD model has the advantage of being able to detect the relationship between driving factors and geographical phenomena without any assumption of linearity. Furthermore, OPGD has four modules, two of which are adopted in this study: the factor detector and interaction detector. We establish the relationship between six factors and positive SAF (LC + SCF + NDVI + SW + SC + ST → positive SAF) and run the OPGD model for each month with a total of 408 runs (34 years × 12 months). Details on the parameter optimization of the OPGD model are provided in the SI. The factor detector uses a q value to quantify the influence of each influencing factor as expressed below: where q is the determinant power of each influencing factor; N and N h indicate the number of samples in the whole region and subregion, respectively; h is the region number (1, 2, 3…); L represents the number of secondary regions; σ 2 is the global variance of Y over the entire study region; σ h is the variance of each subregion h; and the Within sum of squares and the total variance in the whole study area are denoted by SSW and TV, respectively. The q value ranges from 0 to 1, indicating that the influencing factor can explain q × 100% of the explained variable (i.e. positive SAF). The greater q values represent the stronger influence of the factors on the spatial heterogeneity of positive SAF. The interaction detector identifies the interactive effects of two covarying spatial variables on the spatial heterogeneity of positive SAF based on the relative importance of interactions calculated by q values of the factor detector (Song et al 2020). The spatial interaction represents an overlay of two spatial explanatory variables. The weakening, enhancement or independence of two spatial variables can be quantified by the interaction detector. Specifically, there are five interaction types that can be examined using the interaction detector, including nonlinear weakening, univariable weakening, bivariable enhancement, independent, and nonlinear enhancement (Song et al 2020) (table 1). The detailed principles of the interaction detector of OPGD are introduced in the SI.

The spatial pattern and heterogeneity of SAF
The spatial patterns of SAF from 1982 to 2015 are presented in figure 2. From January to March, most Arctic land areas exhibit negative SAFs except for Baffin Island and high latitudes of Greenland in January, YaNAO, parts of Krasnoyarsk Krai and Alaska in February and Chukotka, Banks Island and Victoria Island in March (figures 2(a)-(c)). In April, a positive SAF is found in Alaska, Chukotka and the Northwest Territories, with feedback strengths ranging from −5% K −1 to 0% K −1 (figure 2(d)). From May to June, all Arctic land shows a positive SAF except for the midland of Greenland in May (figures 2(e) and (f)). In particular, high positive SAFs are found in Chukotka, eastern Yakutia, Alaska, Northwest Territories, Banks Island and Victoria Island in June, with feedback strengths above −5% K −1 (figure 2(f)). Notably, in July, most Arctic land areas show a negative SAF, especially in the Krasnoyarsk Krai and Baffin Islands, with a large feedback strength above 10% K −1 (figure 2(g)). From August to December, the positive SAF dominates most regions of Arctic land, with feedback strength ranging from −5% K −1 to 0% K −1 (figures 2(i)-(l)). Overall, a large spatial difference in positive SAF is observed for each month, with CV values beyond 61.5% (figure 2). The largest spatial heterogeneities of positive SAF are found in winter (from December to February), with CV values ranging from 643.9% to 2462.2% (figures 2(a), (b) and (l)). Despite fluctuations, the CV of positive SAF shows an overall decreasing trend at a rate ranging from 0.1% to 3.1% per decade except for February, March, June and November ( figure 3). In addition, the Mann-Kendall test shows that the temporal changes in CV are statistically significant (p < 0.05) for October and December, which decrease at rates of 2.0% and 3.1% per decade, respectively (figure 3). Figure 4 shows the determinant power (q) of six factors on the spatial pattern of positive SAF for each month. Overall, the LC and SCF are the dominant factors determining the spatial pattern of positive SAF in all seasons (figure 4). Specifically, from January to March, LC, SCF and SW have relatively higher determinant power on the spatial pattern of positive SAF, with average q values of 0.22-0.34, 0.14-0.39 and 0.13-0.32, respectively (figures 4(a)-(c)). In addition to LC and SCF, the NDVI also exerts a large influence on the spatial pattern of positive SAF from April to May, with mean q values above 0.32 (figures 4(d) and (e)). During early summer (from June to July), all selected factors have relatively weak determinant power (q ⩽ 0.08) for the spatial pattern of positive SAF (figures 4(f) and (g)). In August, the LC and SCF again have strong contributions to the spatial pattern of positive SAF, with q values of 0.38 and 0.39, respectively ( figure 4(h)). In September, SCF has the highest determinant power, followed by LC, NDVI, SW, ST and SC (figure 4(i)). From October to December, the SCF is the key factor controlling the spatial pattern of positive SAF, with average q values ranging from 0.09 to 0.30 (figures 4(j)-(l)).

Attribution of the spatial heterogeneity of positive SAF
Changing trends of the determinant power (q) for six factors in different months from 1982 to 2015 are summarized in table 2. A decreasing trend is detected in all six factors regarding determinant power in winter (from December to February) (table 2). More    Figure 5 shows 15 pairs of interactions between the six factors affecting the spatial pattern of positive SAF. Specifically, in January, the LC∩SCF and LC∩SW have the largest BE and IND interaction effects, respectively, with an average value of 0.66 ( figure 5(a)). From February to April, the NE and BE interactions between LC and NDVI are the dominant interaction types for the spatial pattern of positive SAF, with average q(LC∩NDVI) values ranging from 0.49 to 0.66 (figures 5(b)-(e)). Similar to the individual effect of a single factor, the interactions of influencing factors have relatively low values in June and July (figures 5(f) and (g)). However, the greatest NE interaction is found between LC and SCF, with q values of 0.17 and 0.19, respectively (figures 5(f) and (g)). In August and September, the BE interaction between SCF and LC has the highest effect on positive SAF, with q(SCF∩LC) values of 0.76 and 0.84, respectively, followed by the interaction of SCF and NDVI, with q values of 0.56 and 0.83, respectively (figures 5(h) and (i)). For October and December, all pairs of interactions have relatively low values compared to other seasons, with the largest NE and UW interactions found for LC∩SW, with q values ranging from 0.26 to 0.28 (figures 5(j)-(l)). More detailed interaction information between different land factors in different months can be found in the SI (tables S1-S12).

Discussion
Changes in snow cover area/extent are recognized as the major drivers of albedo changes in the Arctic (Webb et al 2021). Consequently, the influence of climate perturbation on snow cover is a key process governing Arctic land warming. Generally, snow has a high albedo, and as snow cover decreases, more land surface, which has a much lower albedo, is exposed (Webb et al 2021). In recent decades, most Arctic land areas experienced a significant (p < 0.05) reduction in snow cover during the study period except for Greenland (figure S3), which further contributed to the amplified warming of Arctic land due to the snow-cover-albedo-warming feedback. In addition, the factor detector suggests a high contribution of SCF causes spatial differences in positive SAF at the annual cycle, which is due to the high correlation (p < 0.05) between SCF and albedo at the spatial scale ( figure S4). In addition to snow cover, this study indicates that the vegetation type and coverage also have high contributions to the spatial heterogeneity of Arctic land warming induced by positive SAF. Changes in vegetation coverage inevitably result in changes in surface albedo, which could further affect the regional climate and energy budget of the Earth (Swann et al 2010, Willard et al 2019). For example, in snow-covered areas such as the Arctic, the growth of leaves in springtime and summertime enhances solar absorption due to the reduction in surface albedo (Kelsey et al 2021), which is the reason for the higher contribution of the NDVI to the spatial heterogeneity of positive SAF in this period (figures S5 and 4). In addition, because leaves play a key role in the transpiration process of the hydrological cycle (Seo and Kim 2021), changes in albedo and transpiration rate between different vegetation types induces climate feedbacks. In combination, these processes could result in varying contributions to the spatial pattern of positive SAF. Consequently, the spatial differences in vegetation types over the Arctic further amplifies the spatial heterogeneity of land warming because different types have different albedos and transpiration rates. From figure 2, a higher positive SAF (above −5% K −1 ) is found across western Yakutia and Alaska in June. These regions are generally dominated by dwarf shrubs (figure S2), which have a higher variability in albedo than other LCs, such as herb barrens, moss and wetlands (Loranty et al 2011). In addition, broad-leaved deciduous trees may invade tundra ecosystems more effectively than boreal evergreen trees under climate warming, and because deciduous broad-leaved trees have a higher transpiration rate than needle-leaf evergreen trees, the climate responses may be different (Swann et al 2010).
SW also plays an important role in the spatial pattern of positive SAF because land surface albedo is strongly related to SW. Previous studies have reported that surface albedo has a linear relationship with soil moisture (Roxy et al 2010, Li et al 2019, Yang et al  2020). Generally, dry soil has a higher albedo (above 0.7), while wet soil has a lower albedo (below 0.5) (Atlaskina et al 2015). As indicated by figure S6, the SW in the Northwest Territories, Nunavut, and Baffin Island is relatively low. Hence, these regions are more likely to be affected by the different types of vegetation, resulting in greater SAF, which is confirmed in figure 2. In addition, the spatial differences in SW also have a great impact on vegetation growth and soil respiration. In general, vegetation growth is highly sensitive to soil moisture, and therefore spatial differences in SW could affect the distribution, composition and abundance of vegetation ( This study quantifies the strengths of seasonal SAF and explores the contributions of different land factors to spatial patterns of positive SAF over Arctic land. Although it is out of the scope of this paper to consider other influencing factors, such as aerosols, solar irradiance, greenhouse gases, sea ice, and clouds, it is worth giving attention to the nonlinear interactions between different land and atmospheric components of the climate system. Indeed, sea ice-albedo feedback, cloud feedback, and water vapor feedback could also amplify or restrain climate perturbations. For example, the feedbacks related to sea ice reduction are shown to be critical to polar amplification (Dai 2021, You et al 2021, which could further accelerate the melting of snow and/or the expansion of vegetation (Mudryk et al 2017). Additionally, the distribution and frequency of rainfall and snowfall are impacted by aerosols, water vapor and cloud feedbacks (Pan et al 2020). More importantly, the deposition of aerosols (e.g. black carbon, dust, and nitrate) also contributes to the reduction of albedo in polluted regions, although the forcing is less compared to the preindustrial period (Alessandri et al 2021). Nevertheless, the deposition of nitrogen and phosphorous may enhance biogeochemical cycles and improve the productivity of terrestrial ecosystems (Mahowald et al 2017). Under this circumstance, the expansion or reduction of vegetation cover could cause feedback by regulating the emission of bioaerosols and their precursors, whereas the release of mineral dust could also be affected because it is preferentially emitted from dry and unvegetated soils (Ganzeveld et al 2010). Moreover, the multiple biophysical interactions among snow, vegetation and soil tend to have larger effects than the absorption of solar radiation. It has been shown that when snow is present, the expansion of forested areas will lead to positive surface net longwave radiation due to increasing longwave radiation, which could cause feedback to seasonal warming, especially in temperate and subarctic areas (Todt et al 2019). The future increase in regional temperature is highly affected by these complex interactions and cumulative effects at various space and time scales.
Representing the dynamics of snow, vegetation and soil, Earth system models (ESMs) are valuable tools for understanding the role of SAF in amplifying or reducing anthropogenic climate warming at the regional scale (Winton 2006). However, a range of results are found among state-of-the-art ESMs due to uncertainties in setting land surface parameter schemes (Thackeray and Fletcher 2016, Andry et al 2017, Thackeray and Hall 2019. As a result, current ESMs must be improved and better constrained to achieve reliable projections of future climate change, especially for warming amplification induced by SAF. Thus, this work could provide a benchmark for improving the representation of the processes governing the variations in seasonal SAF in the next generation of ESMs. More importantly, the impacts of different land parameters and their interactions have large implications for optimizing relevant parameterizations of ESMs.

Conclusions
In this study, the strengths of seasonal SAF are quantified for the period between 1982 and 2015 in the Arctic. A large positive SAF (above −5% K −1 ) is observed across Baffin Island in January and eastern Yakutia in June. In contrast, a large negative SAF (beyond 5% K −1 ) is estimated in Nunavut, Banks Island, Victoria Island, Chukotka and low latitudes of Greenland in January and Nunavut, Baffin Island and Krasnoyarsk Krai in July. As indicated by the CV, large spatial differences in positive SAF are found in the Arctic, which contribute to the high spatial heterogeneity of Arctic land warming. Furthermore, based on the OPGD model, the determinant power and the interactive effects of six factors are examined. Generally, the rank of determinant power is SCF (0.26 ± 0.14) > LC (0.23 ± 0.13) > NDVI (0.18 ± 0.12) > SW (0.13 ± 0.08) > SC (0.07 ± 0.03) > ST (0.06 ± 0.03).
In addition, the interactions between SCF, LC and SW exhibit a further influence on the spatial heterogeneity of Arctic land warming induced by positive SAF in March, June and July.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files). Data will be available from 02 August 2021.