Underestimation of the impact of land cover change on the biophysical environment of the Arctic and boreal region of North America

The Arctic and Boreal Region (ABR) is subject to extensive land cover change (LCC) due to elements such as wildfire, permafrost thaw, and shrubification. The natural and anthropogenic ecosystem transitions (i.e. LCC) alter key ecosystem characteristics including land surface temperature (LST), albedo, and evapotranspiration (ET). These biophysical variables are important in controlling surface energy balance, water exchange, and carbon uptake which are important factors influencing the warming trend over the ABR. However, to what extent these variables are sensitive to various LCC in heterogeneous systems such as ABR is still an open question. In this study, we use a novel data-driven approach based on high-resolution land cover data (2003 and 2013) over four million km2 to estimate the impact of multiple types of ecosystem transitions on LST, albedo, and ET. We also disentangle the contribution of LCC vs. natural variability of the system in changes in biophysical variables. Our results indicate that from 2003 to 2013 about 46% (∼2 million km2) of the region experienced LCC, which drove measurable changes to the biophysical environment across ABR over the study period. In almost half of the cases, LCC imposes a change in biophysical variables against the natural variability of the system. For example, in ∼35% of cases, natural variability led to −1.4 ± 0.9 K annual LST reduction, while LCC resulted in a 0.9 ± 0.6 K LST increase, which dampened the decrease in LST due to natural variability. In some cases, the impact of LCC was strong enough to reverse the sign of the overall change. Our results further demonstrate the contrasting sensitivity of biophysical variables to specific LCC. For instance, conversion of sparsely vegetated land to a shrub (i.e. shrubification) significantly decreased annual LST (−2.2 ± 0.1 K); whereas sparsely vegetated land to bare ground increased annual LST (1.6 ± 0.06 K). We additionally highlight the interplay between albedo and ET in driving changes in annual and seasonal LST. Whether our findings are generalizable to the spatial and temporal domain outside of our data used here is unknown, but merits future research due to the importance of the interactions between LCC and biophysical variables.


Introduction
The warming rate of the Arctic and Boreal Region (ABR) is twice that of the rest of the world and will likely remain this way for the rest of the 21st century independent of the shared socio-economic pathways (Stocker et al 2013, Masson-Delmotte et al 2021).This phenomenon is known as Arctic amplification (Previdi et al 2021) and contributes to the increasing intensity of global warming.One of the main consequences of Arctic amplification is intensified land cover change (LCC), such as shrub expansion, loss of forest, changes in the dominant forest types, and various other LCC types mainly caused by climate change and disturbance (Tape et al 2006, Soja et al 2007, Carpino et al 2018, 2021, Wang et al 2020, 2021).Biophysical variables such as land surface temperature (LST), albedo, and evapotranspiration (ET) are among the major components of surface energy balance and are highly sensitive to LCC (Lean and Warrilow 1989, Snyder et al 2004, Swann et al 2010, Alkama and Cescatti 2016, Liao et al 2018).LCC alters these variables locally, which ultimately results in feedback to the global climate system that in turn influences the structure and function of the ABR.A better understanding of the interaction between LCC and biophysical variables is of critical importance to ecological and climate change forecasts (Feddema et al 2005, Bala et al 2007), carbon sink and source studies, natural climate solution implementation (Field andMach 2017, Griscom et al 2017), and local and regional land management decision making (Duveiller et al 2020).
The ABR is an ideal study region to quantify the impacts of LCC on biophysical variables due to its historic above-average warming and multiple known ecosystem transitions such as the expansion of deciduous forests (DFs), reduction in evergreen forests (EFs) (Bonan et al 1992, Swann et al 2010), shrubification (Mekonnen et al 2021), shrinking wetlands in one area, and expanding in another area (Kreplin et al 2021).Despite various active LCC in the region, most studies on the sensitivity of biophysical variables to LCCs are limited to some major ecosystem transitions such as deforestation and afforestation (Swann et al 2010, Lee et al 2011, Zhao and Jackson 2014, Li et al 2015, Alkama and Cescatti 2016, Chen and Dirmeyer 2020, Novick and Katul 2020) or a few specific ecosystem transitions (Claussen et al 2001, Bathiany et al 2010, Dass et al 2013, Devaraju et al 2015).However, there are many types of LCC, which can impact biophysical variables LST, albedo, and ET differentially, leading to both positive and negative feedbacks to the climate system.Due to limited data availability and the relative inaccessibility of the region, there has yet to be a full evaluation of multiple interacting ecosystem transitions emerging across the region that extends over space and time.
LCC impacts LST, albedo, and ET through a series of complex processes such as light interaction with canopy and understory, surface and aerodynamic resistance, photosynthesis, and other ecophysiological processes (Bonan 2019, Chen and Dirmeyer 2020, Zhang et al 2020).The extent of influence of each process mainly depends on the type of ecosystem transition and the local background climate (Li et al 2016, Forzieri et al 2017, Novick and Katul 2020).For example, a transition from forest to shrub in cold regions has been found to result in reduced LST, since cooling due to increased albedo is generally dominant over warming due to decreased ET (Bonan et al 1995).In contrast, this same transition from forest to shrub in tropical regions has been found to result in increased LST, since warming associated with decreased ET is the dominant factor (Li et al 2015).Better quantification of potentially contrasting effects of different types of ecosystem transitions on LST, albedo, and ET is needed to more fully predict the net impact of a given ecosystem transition.
The most common approach to estimate the impact of LCC on biophysical variables is based on attribution models (Juang et al 2007, Rigden and Li 2017, Liao et al 2018, Chen and Dirmeyer 2020).These are analytical models based on surface energy balance and are mostly applied to LST but can be extended to other variables such as ET (Wang et al 2019a).While these models are powerful and can physically attribute changes in biophysical variables to the governing processes, they require detailed information such as surface and aerodynamic resistance that are usually provided by eddy covariance (EC) towers or Earth system models (ESMs).However, in ABR the distribution of EC towers is sparse due to remoteness and challenges with instrumentation (Pallandt et al 2021).The application of ESMs is also limited due to a variety of reasons including the lack of proper representation of land cover in models (Fisher and Koven 2020) and poor partitioning of available energy into latent and sensible heat fluxes (Pitman et al 2009, de Noblet-ducoudr et al 2012).Because of these limitations, ESMs have been found to vary widely in their predictions of LCC impact on the biophysical environment (Duveiller et al 2018).For example, Pitman et al (2009) analyzed seven different ESMs and showed that the difference in temperature change due to LCC between models can be up to 4 • C.
Remote sensing is a powerful technique to overcome previous limitations and estimate the impact of LCC on LST, albedo, and ET.There are multiple studies on the interaction between LCC, mostly forest ecosystems, and biophysical variables using satellite observations (Peng et al 2014, Zhao and Jackson 2014, Li et al 2015, Forzieri et al 2017, Duveiller et al 2018).The work done by Alkama and Cescatti (2016), is one of the well-recognized studies on the sensitivity of biophysical variables to afforestation and deforestation.Their methodology disentangles the changes in LST due to natural variability vs. changes imposed by deforestation and afforestation.While their work is limited to forest ecosystems, it sets a robust framework that can be extended to other land cover transitions.To our knowledge, Duveiller et al (2018) provided one of the earliest works on the sensitivity of biophysical variables to multiple types of ecosystem transitions.They used a static land cover map with 300 m spatial resolution and the space-for-time method to calculate the potential change in surface energy balance components when an LCC happens.The coarse resolution of the land cover map used in Duveiller et al (2018) imposes some challenges in interpreting their results over ABR as the heterogeneity of this region leads to many land cover transitions over small areas.Landsat 30 m resolution provides relatively high spatial resolution data that recently have been used to map detailed land cover maps over the region (Wang et al 2020).This dataset along with widely validated biophysical products from the Moderate Resolution Imaging Spectroradiometer (MODIS) can be used to get more insight into interactions between LCC and biophysical variables.
Here we advance previous efforts by integrating land cover data with 30 m spatial resolution over approximately four million km 2 capable of resolving fine-scale LCC, and LST, albedo, and ET observations with 1000 m spatial resolutions to quantify the impact of the full range ecosystem transitions occurring across ABR.The main question of this study is what is the impact of different types of ecosystem transitions on albedo, LST, and ET over the ABR?In our study, we disentangle the changes in biophysical variables due to natural variability and LCC.We build upon previous works and consider more detailed ecosystem transitions.

Study area and datasets
The study region is the core domain of the NASA Arctic-Boreal Vulnerability Experiment project which includes most of Alaska and a large portion of western Canada (figure 1) and is slightly more than four million km 2 .In this study, we refer to this domain as ABR.The study period includes the years 2003 and 2013.The main reason for choosing this period is the overlap between different data used in this study.The datasets include a Landsat-derived land cover map (Wang et al 2019b), MODIS-based blue sky albedo (Solvik et al 2019), and the global actual ET dataset (Zhang et al 2019), and LST (Hulley 2017).Below we describe the preprocessing involved for each product.The analyses are mostly done using a series of open-source packages developed for the Python language in addition to other open-source software such as Google Earth Engine (Gorelick et al 2017) and QGIS.The code and datasets to reproduce the results of this study can be found in: https:// github.com/hamiddashti/lcc-biophys-arctic.

Data preparation 2.2.1. Land cover data
The land cover map used in this study was recently developed for the ABR using Landsat (Wang et al 2020).The spatial resolution of the land cover data is 30 m and it covers the 1984-2014 period.The data can be downloaded through the ORNL DAAC (10.3334/ORNLDAAC/1691).This product has ten different land cover types (table S1) including EF, DF, shrublands (shrub), herbaceous (herb), sparsely vegetated (sparse), barren, fen, bog, shallow/littoral, and water classes (Wang et al 2019b).The dataset is in Canada Albers Equal Area Conic (EPSG: 102 001) projection system.

LST
In this study, we used the LST product, version 6 (MYD21v006) estimated from the MODIS instrument onboard the Aqua satellite.The data were downloaded from the Application for Extracting and Exploring Analysis Ready Samples (AρρEEARS, https://lpdaac.usgs.gov/tools/appeears)provided by the US Geological Survey.The data was ordered in Albers equal-area coordinate system.AρρEEARS uses the nearest neighbor resampling method when converting between different coordinate systems.The delivered data has a pixel size of ∼1000 m and an eight days temporal resolution.In the next step, we filtered data based on the quality flags.Pixels with LST accuracy flagged as 'good performance' (LST accuracy 1-1.5 K) and 'excellent performance' (LST accuracy <1 k) are kept in the analyses.MODIS daytime and nighttime LST were recorded at approximately 13:30 and 1:30 local time, respectively.This time is close to maximum and minimum temperature (Duveiller et al 2018).We calculated the mean LST by taking the average between day and night estimates.Finally, the eight days data were resampled to annual by taking the average of all measurements in each year.

Blue sky albedo
The albedo data are based on the recently developed MODIS-based daily albedo with 500 m pixel size for Alaska and Canada (Solvik et al 2019).Most available albedo datasets are modeled albedo under pure diffuse illumination (white sky albedo) or pure direct illumination (black sky) (Potter et al 2020).The dataset used in this study is 'blue sky' albedo which refers to the albedo under actual sky illumination conditions and is appropriate for use in higher latitudes with large solar zenith angles covered by highly reflective materials such as snow and ice.Our quality control is the same as (Potter et al 2020).In the next step, we take the annual mean of albedo values and then reproject and resize data into LST pixel size using a bilinear resampling method.Note that over higher latitudes elements such as polar nights impact annual albedo.However, we suspect that polar night will not have a dramatic impact on our results, since regardless of land cover type, during this time there is no reflected radiation, and surface albedo no longer matters.The polar night is also relatively short (about a month), giving us multiple months of winter with good albedo data which contribute to the annual albedo.

ET
The Penman-Monteith-Leuning (PML-V2) ET dataset is used in this study (Zhang et al 2016(Zhang et al , 2019, Gan et al 2018).The algorithm used in the PML_V2 is based on biophysical processes such as canopy conductance and carbon constraint on ET.It takes MODIS data (leaf area index, albedo, and emissivity) together with Global Land Data Assimilation System (GLDAS) meteorological forcing data as model inputs.The spatial resolution of the PML-V2 ET product is 500 m and its temporal resolution is eight days.PML-V2 is well-calibrated against eight days measurements at 95 widely distributed flux towers (Zhang et al 2019).PML-V2 product provides different components of ET including evaporation from the soil, transpiration from the plant canopy, evaporation of precipitation intercepted by the vegetation, and evaporation from water bodies, snow, and ice.The code of the PML-V2 model is available at https://github.com/gee-hydro/gee_PML.We extracted monthly and annual data of all ET components in the study area for years 2003 and 2013 by summing up all eight days products in each month and year using Google Earth Engine.The processing code is publicly available at https://code.earthengine.google.com/17ded6f2e121b942f54d1ed75ef9f050.

Disentangling natural variability from LCC impact
The annual and seasonal mean values of biophysical variables (section 2.2) were used.In our analyses, we only included those pixels that we have high-quality LST, albedo, and ET data.Our annual data might be biased towards the growing season as we have fewer wintertime high-quality data, in particular for albedo.However, we refer to them to these data as annual since we still have some wintertime data.The method used in the study to disentangle the impact of natural variability from LCC on biophysical variables is based on Alkama and Cescatti (2016).We assume the change in a given biophysical variable from 2003 to 2013 has two components: (a) the change caused by LCC, and (b) the change due to the natural variability of the system.We can represent these two components in equation ( 1) below: where X refers to a given biophysical variable (i.e.LST, albedo, and ET).∆X Total is the total change in a biophysical variable from 2003 to 2013, ∆X LCC is the changes in biophysical variables caused by LCC, and ∆X NV is the natural variability of a biophysical variable.In this study, we are most interested in ∆X LCC , which can be calculated as: From observations, we already have where d k is the distance of the kth unchanged pixel to the pixel of interest and the ∆X Total, K is the total difference in the biophysical variable from 2003 to 2013 for the pixel.(e) Finally, we calculate the ∆X LCC of the central pixel by subtracting the associated ∆X NV from ∆X Total of the pixel.(f) Move to the next 1000 m pixel for a given biophysical variable and repeat the process if the given pixel experienced an LCC of more than 2%.

Impact of LCC on biophysical variables
To calculate the impact of LCC on biophysical variables we used a regression analysis developed by Alkama and Cescatti (2016).To explain the approach, let us assume we are interested in the impact of land cover type 1 (LC1) to land cover type 2 (LC2) on ∆X LCC .regressions between LCC bins and ∆X LCC using bootstrapping with replacement with the weights representing the number of observations in each bin.We extract the mean and standard deviation of the estimated regression parameters as the final estimates and their uncertainty, respectively.If the extent of LCC from or to a given land cover type was less than 50%, we consider that regression unreliable and excluded it from analyses.This was a common situation with ecosystem transitions to or from water, littoral, and bogs, and thus we did not consider them in our final analyses.We note, that more detail on this method can be found in Alkama and Cescatti (2016).
The regression models have the following forms: where β 0 and β 1 are the intercept and slope.The intercept can be interpreted as the expected change of ∆X LCC if LCC LC1−> LC2 is zero.The intercept value should be zero, since if there is no change in land cover then there is no change in a biophysical variable due to LCC.However, in practice, this value is not exactly zero since neither the data nor our approach is perfect.A small value of intercept is expected.The slope parameter is the expected change in ∆X LCC for a unit change of LC1 to LC2.Thus, the slope values can be considered as the sensitivity of ∆X LCC to LCC which can be read as the impact of LCC on ∆X LCC .The conceptual figure S1 depicts the steps explained above.In our study, we disentangle the contribution of LCC vs. natural variability in total changes of LST, albedo, and ET.Overall, the region experienced a decrease in LST (−1.1 ± 1.7 K), a slight decrease in albedo (−0.008 ± 0.028), and an increase in ET (9.2 ± 27.3 mm) (table 2) from 2003 to 2013.The resulted decrease in LST is a function of interannual variability and we did not observe any significant cooling trend from 2003 to 2013.∆X NV has a more regional and uniform pattern compared to the ∆X LCC (figure S1).In all cases, the regional mean of ∆X NV signal is larger than the regional mean of ∆X LCC .For example, the change in LST due to the natural variability of the system and LCC are −1.0 ± 1.2 and −0.1 ± 1.13 K, respectively.Moreover, the range of change in biophysical variables due to LCC is more symmetrical than natural variability.About 8% of the pixels with a change in LST due to LCC experienced a more extreme increase in LST (∆LST LCC ⩾ 1.5 K) and ∼11% of pixels experienced an extreme decrease in LST (∆LST LCC ⩽ −1.5 K) (figure S2).

Natural variability vs. LCC
Disentangling ∆X LCC and ∆X NV based on the direction of change shows that the magnitude of the effect of LCC on biophysical variables can be larger than the natural variability and can dramatically alter the sign or magnitude of increase or decrease (figure 2).In approximately half of the areas that experienced LCC, change in biophysical variables due to LCC has an opposite sign compared to change due to the natural variability.About 45% of the changed area experienced the opposite signs in ∆LST LCC and ∆LST NV and ∼53% and ∼49% for albedo and ET, respectively.In 34.5% of the area, LCCs imposed an increase in LST, which dampened the observed decrease in ∆LST NV by about 0.9 ± 0.6 K.In contrast, in 43.5% of the changing area, an LCC intensified the ∆LST NV by about 1.0 ± 0.7 K. Interestingly, when considering changes in albedo, we found that the reduction in albedo due to natural variability was completely ameliorated by the increase in albedo driven by LCC, resulting in a change of ∆albedo Total sign.A similar pattern can be seen for LST where the ∆LST NV increased by 0.6 ± 0.4 K but the ∆LST LCC decreased by −0.9 ± 0.6 K leading to a −0.2 ± 0.8 K decrease in ∆LST Total (figure 2).

Impacts of major ecosystem transition
We estimated the impact of ecosystem transitions observed over the region on biophysical variables.We focused on ecosystem transitions where the extent of LCC is more than 50% to highlight the dominant transitions where we have more confidence in their estimated slope.Figure 3 shows the impact of total gain or loss of a land cover, regardless of the type of the transition, on ∆X LCC and figure 4 shows the impact of specific LCC on ∆X LCC for those ecosystem transitions that we have a reliable estimate.Specific types of ecosystem transitions were found to lead to different responses in biophysical variables and sometimes these responses have opposite signs.For In each panel, the red arrows show the increase or decrease of these variables.The area is the fraction of pixels that experienced a decrease and increase in ∆X components associated with each subplot.We can observe that for all three variables in almost half of the cases the ∆XLCC has a sign opposite to the ∆XNV.Taking the mean of all these changes (as in table 2) leads to an underestimation of the impact of LCC on the biophysical environment of the region.
instance, the conversion of EF to sparse vegetation led to a decrease in both albedos (−0.009 ± 0.001) and ET (−32.7 ± 0.1 mm), which increases LST (+0.9 ± 0.0 K).Conversely, the transition of EF to shrub led to an increase in albedo (0.025 ± 0.004) and a decrease in ET (−45.4 ± 1.0 mm), which resulted in a decrease in LST (−1.0 ± 0.0 K).This contrasting behavior is noticeable as one expects the conversion of forest to sparse or shrub should lead to the same albedo response.However, we observe in practice the sensitivity of forest albedo to sparse transition is less than transition to shrub.
Shrub gain from 2003 to 2013 was found to result in a decrease in LST, an increase in albedo, and a decrease in ET (figure 3).Focusing on the individual ecosystem shrub-sparse, shrubfen, and shrub-herb, we observe a more interesting pattern.The sensitivity of the biophysical variables to the shrub-herbaceous transitions is negligible (figure 4).The transition of sparse to shrub and fen to shrub led to a decrease in LST (−2.2 ± 0.1 and −1.0 ± 0.1 K, respectively) and an increase in albedo (0.058 ± 0.006 and 0.036 ± 0.007, respectively).However, the response of ET for these transitions is in the opposite direction; an increase in ET (38.3 ± 2.9 mm) for sparse to shrub and a decrease in ET (−10.7 ± 0.5 mm) for fen to shrub.
The fen gain led to a slight reduction in ET.Considering fens are wetlands, this reduction might look surprising.However, a closer investigation shows that most transitions of ecosystems to fen led to an increase in ET, except for deciduous and EF transitions (figure S4).Thus, considering that forest ET is generally higher than fen over ABR (table 1), the reduction of ET following fen gain is expected.In the transition of a shrub to fen, ET increased (10.6 ± 0.52 mm), but albedo decreased (−0.036 ± 0.007), possibly due to lower water albedo, and ultimately LST increased by 1 K. Overall, focusing on total gain and loss for different land cover types, Black lines are the mean of 10 4 bootstrap regression models.The estimated slopes reported are the mean of bootstrap regressions ± standard deviations and can be interpreted as the impact of a land cover gain and loss on a biophysical variable.The intercepts are very close to zero, hence not shown.The gain and loss analyses can be interpreted as the mean impact of all types of ecosystem transitions on the biophysical environment.
We also investigate the seasonal sensitivity of ∆LST LCC to various LCCs (figures 5 and S4).Due to the lack of enough wintertime albedo for the ABR, we did not perform these analyses for the albedo and ET.The difference in the seasonal impact of LCC on ∆LST can be observed.Notably, in all major ecosystem transitions, the direction of change in LST during winter and summer has an inverse sign.The transition of EF to all classes except DF leads to wintertime cooling and summertime warming.Spring follows the same decrease and increases pattern as wintertime however it is less intense.For example, Loss of EF and DF on average leads to a ∼−1.5 K decrease in LST during winter while it is ∼−1.0K during springtime.In most transitions, fall follows summertime with a few exceptions such as transitions of shrubs to herbaceous and fen that leads to ∼0.5 and

Discussion
Our analyses indicate that from 2003 to 2013 approximately 45% of the region (∼2 million km 2 ) experienced some type of LCC which has already caused changes to the biophysical environment across the Arctic Boreal zone.It is common in the literature to report one value for the change in biophysical variables due to LCC for this region (Alkama and Cescatti 2016, Forzieri et al 2017, Perugini et al 2017, Duveiller et al 2018, 2020).However, our findings indicate that reporting the mean without consideration of the sign of change relative to natural variability may lead to the underestimation of the influence of LCC on biophysical variables.For example, from table 2 one can conclude that the −0.1 K ∆LST LCC is ten times smaller than −1.0 K ∆LST NV over the entire region.In another word, LCC played a negligible role in changing LST from 2003 to 2013.However, when considering the sign of the change (figure 2), we can see the change in ∆LST LCC is as large as ∆LST NV .Indeed in ∼11% of the cases the sign of ∆LST Total changed from positive to negative due to the dominant impact of the ∆LST LCC compared to ∆LST NV .Thus, the influence of LCC is significant and can either intensify or lessen changes in the biophysical environment due to natural variability.Taking the mean over the entire region leads to the underestimation of LCC's impact on biophysical variables due to the cancelation of these opposing signs.Whether ecosystems are negatively or positively impacted by LCCdriven negation or intensification of the natural variability of the biophysical environment is contextual and depends on the location.
Our findings highlight the importance of carefully considering multiple land cover transitions.The contrasting response of biophysical variables to different types of LCC is of particular interest.For instance, several studies show that deforestation over the Arctic and Boreal should increase albedo and hence decrease LST (Randerson et al 2006, Li et al 2015).Our results are consistent with these studies, especially for DF where deforestation decreased LST (figure 3).However, when considering the specific types of deforestation, we get a contrasting response.For example, EF transition to the sparsely vegetated decreased albedo and increased LST (figure 4).By contrast, EF to shrub increased albedo and decreased LST.In our land cover data, sparse is defined as pixels with '10%-30% canopy coverage and rock underneath' (table S1).Thus, in EF to sparse transition, the summertime darker rocks or soil may contribute more to the reduction of albedo compared to early spring/late fall increase in albedo due to snow.And in the case of EF to shrub, shrubs may cover the lower albedo soil and rock hence less profound impact of summertime decrease in albedo.When aggregating these two transitions as 'deforestation' the stronger increase of albedo in EF to shrub (0.025 ± 0.004) might dominate the small decrease of albedo in EF to sparse (−0.009 ± 0.001) leading to an increase in overall albedo.Another interesting observation is the greater reduction of ET in EF to shrub than EF to sparse conversion from 2003 to 2013.It is not truly clear to us what exactly causes the ET responses to these ecosystem transitions.However, the most plausible explanation can be attributed to different components of ET.For example, while it is true sparse vegetation has smaller transpiration compared to shrubs due to less vegetation, the soil evaporation might compensate and probably dominate the total ET.Indeed, there is evidence that moss which is common in sparsely vegetated areas can evaporate more than in open waters (Heijmans et al 2004).The conversion of forest-to-sparsely vegetated and potential warming is more common than some cooling transitions such as forest-to-shrub, at least over the decade that we considered here (table S2).Thus, to better understand the net impact of forest land cover conversion including reforestation and afforestation events the type of transition and its resulting biophysical impact must be accounted for.
Notably, we found that shrub gain (i.e.shrubification) has generally driven a decrease in LST.Annual albedo and ET were found to increase due to shrub gain, which contributed to an overall annual decrease in LST (figures 3 and 4).However, it becomes clear from our seasonal analysis of LST that there are contrasting seasonal effects of shrubification (figure 5).Considering the most common shrubification pathway, sparse-to-shrub (Cornelissen et al 2001, Elmendorf et al 2012, Hollister et al 2015, Mekonnen et al 2021), we found strong summer (−3.7 ± 0.2 K) and fall (−1.0 ± 0.1 K) cooling that was dominant over winter (1.1 ± 0.0 K) and spring (0.3 ± 0.0 K) warming.While a seasonal analysis was not possible for albedo and ET, we speculate that increases in albedo and ET during summer and fall were the major driving factors behind annual decreases in LST due to shrubification from 2003 to 2013.Increases in albedo due to sparse-to-shrub transition are expected when soils are dark; whereas an increase in ET due to sparse-to-shrub transition is expected in all cases.Thus, the summer/spring decrease over the entire region due to shrubification might be more attributable to the increase in ET.
At first glance, our findings on shrubification might be in contradiction with previous studies that have reported a positive feedback between shrubification and climate warming (Bonfils et al 2012, Loranty and Goetz 2012, Mekonnen et al 2021).However, we should note that some fundamental differences that make the direct comparison challenging.Most studies linking shrubification to positive feedback on the warming are trend analyses based on long-term data (Beringer et al 2005, Chapin et al 2005, Loranty and Goetz 2012, Lafleur and Humphreys 2018) or ESMs (Chapin et al 2000, Bonfils et al 2012).Our proposed method is not designed for trend analyses, rather it decomposes changes in biophysical variables between the natural variability of the system and the impact of various LCC on these variables between two end years at the regional scale.Further, the word warming in most trend analysis studies refers to the air temperature rather than LST used in this study.The relationship between LST and air temperature is nonlinear and complex (Novick andKatul 2020, Zhang et al 2020).To what extent the changes in LST due to LCC lead to changes in air temperature remains an open question.And finally, the common assumption in analyses based on ESMs is that albedo is the dominant factor controlling LST (Juang et al 2007, Swann et al 2010).However, we showed that it is likely that changes in ecophysiological processes controlling ET dominate the albedo's impact on LST, at least over the study period.These fundamental differences in variables used, method and assumptions are a common occurrence in studies concerning Arctic greening and its implications (Myers-Smith et al 2020).Thus, we encourage careful consideration when comparing our results with other studies.
It has been shown that the temperature seasonality-conventionally defined as the difference between summer and winter temperatures-has been diminished over the Arctic (Boreal) region equivalent to 4 • to 7 • latitudinal shift equatorward (Xu et al 2013).This phenomenon has major ecological, climatic and societal impacts.We showed that LST change has the opposite sign in all major ecosystem transitions (figure 5), hence intensifying the summerwinter temperature differences.Thus, LCC might be a pathway that can help to reduce the diminishing temperature seasonality over the region.The seasonal response of biophysical variables can also shed light on some of the unexpected LST observations reported in this study.As an example, we discussed the possible dominant role of summer and fall ET in driving decreases in LST in the case of the sparse-to-shrub ecosystem transition.In contrast, our results imply that summertime albedo might be a dominant factor in driving increases in LST in the shrub-to-fen transition (figure 5).Depending on location and other environmental factors, fens during winter might be covered by ice and snow or covered with water.Both covers (snow/ice and water) have sharp and opposite albedo responses.During summer, the transition of a shrub to fen leads to an increase in ET but a decrease in albedo due to the presence of water.Ultimately, these complex seasonal interactions led to an overall reduction in albedo with an increase in LST during winter, spring, and summer, and these seasonal changes appear to drive the overall increase in ET with a decrease in LST.Climate change and Arctic amplification are rapidly changing ice and snow cover (Liu et 2020), and other environmental factors.There is little understanding of the seasonal impacts of these changes on biophysical variables and ultimately surface energy balance over the region.
Our work provides a new perspective on using high to moderate resolution data to study the impact of LCC on the biophysical environment.However, there are some limitations to this approach that we would like to highlight to encourage future research aimed at methodological advances.First, a different selection of the two end years might lead to slightly different results.As we discussed in the introduction the impact of LCC on biophysical variables depends on the background climate.As an example, the conversion of forest to sparse in a dry year vs. wet year may lead to a different response of ∆LST LCC .However, this is an inherent characteristic of the system and not necessarily a flaw of the current study.The generalization of results presented here outside of the spatial and temporal domain of this study remains an open question.Second, in this study, we used data from the entire ABR domain, and consequently, our results need to be interpreted at the regional scale.In another word, the reported sensitivities are the mean of changes in biophysical variables due to various LCC over the entire ABR, rather than local sites.It is important to note that local factors such as latitudinal differences, species, diseases, and topography might impact the response of biophysical variables to a given LCC.It was not practical to consider all these elements since disentangling the region further based on biogeography, latitude, and other factors would result in a reduction of data points and statistical power.However, we were able to attribute some ∆LST LCC extremes to local factors such as wildfire (figure S2).We found that in many cases in central Alaska and eastern ARB, there is an overlap between ∆LST LCC extremes and wildfire scars (generated by Loboda et al 2017) from 2003 and 2013.A promising pathway for future studies is to incorporate higher resolution data to resolve local factors mediating the relationship between LCC and the biophysical environment.Finally, in this study, we only considered LST, albedo, and ET and not the full components of surface energy balance.It would be interesting to investigate the impact of LCC on other components, especially during wintertime when factors such as polar nights play important role in surface energy budget.
To limit global warming to below the 2 • C target set by the Paris Agreement (UNFCCC 2015), improved ecosystem stewardship, or the responsible use and protection of the natural environment, has been identified as an important mitigation approach (Field and Mach 2017, Griscom et al 2017, Fargione et al 2018).Key pathways toward improved ecosystem stewardship include natural climate solutions that take into account the impact of LCC, such as reforestation, afforestation, grassland conservation, and fire management (Griscom et al 2017, Fargione et al 2018).Here, we show that biophysical variables have different and sometimes contrasting responses to ecosystem transitions which might lead to the underestimation of the impact of LCC on the biophysical environment.We further put forward some baseline estimates of LCC impacts on key biophysical variables LST, albedo, and ET for the rapidly changing Arctic Boreal region, and present an approach that could be replicated for other regions globally.

Conclusion
The LST, albedo, and ET play an important role in regulating the energy balance in the Arctic.In this study, we moved beyond the general consideration of common LCC pathways such as deforestation and shrubification and quantify the sensitivity of these key biophysical variables to a wide range of ecosystem transitions.We showed that in half of the cases the LCC imposes a change in the LST, albedo, and ET in an opposite direction of the natural variability of the system.Thus, precaution is needed when aggregating these changes in biophysical variables over large regions such as the ABR, since the opposite signs of natural and LCC impacts will result in cancellation and an overall underestimation of the impact of LCC.Moreover, we also clearly demonstrate that the type of LCC transition matters.For example, deforestation and afforestation are not enough to explain the complex interaction of forest change with biophysical variables.Conversion of forest to shrub and sparsely vegetated lands are both deforestation but the former decreases and the latter increases LST.Finally, we highlight some of the limitations of this study such as the need for higher resolution data, considering only end years (2003 and 2013) rather than all years, and not incorporating the full surface energy balance components.This study is a first step toward a deeper understanding of the relationship between LCC and the biophysical environment.

Figure 1 .
Figure 1.Land cover and land cover change (LCC) across the Arctic and Boreal Region (ABR) between 2003 and 2013.The dominant land cover in (A) 2003 and (B) 2013; (C) the maximum fractional change of LCC that happened between the two end years; and (D) an example of various LCC over a relatively small area possibly due to succession following a disturbance such as wildfires.

Figure 2 .
Figure 2. The mean of (panel A) ∆LST (K), (panel B) ∆Albedo, and (panel C) ∆ET (mm) for the region with different increased and decreased directions in ∆XLCC and ∆XNV.In each panel, the red arrows show the increase or decrease of these variables.The area is the fraction of pixels that experienced a decrease and increase in ∆X components associated with each subplot.We can observe that for all three variables in almost half of the cases the ∆XLCC has a sign opposite to the ∆XNV.Taking the mean of all these changes (as in table 2) leads to an underestimation of the impact of LCC on the biophysical environment of the region.

Figure 3 .
Figure 3.The impact of total gain and loss of different land cover types on (panel A) LST, (panel B) albedo, and (panel C) ET.Black lines are the mean of 10 4 bootstrap regression models.The estimated slopes reported are the mean of bootstrap regressions ± standard deviations and can be interpreted as the impact of a land cover gain and loss on a biophysical variable.The intercepts are very close to zero, hence not shown.The gain and loss analyses can be interpreted as the mean impact of all types of ecosystem transitions on the biophysical environment.

Figure 4 .
Figure 4.The estimated slopes are associated with the full range of measurable ecosystem transitions for (A) land surface temperature (LST), (B) albedo, and (C) evapotranspiration (ET).The values shown are for ecosystem transitions from the year 2003 (y-axis) to the year 2013 (x-axis).Since relationships are linear, the transitions associated with the lower triangle has a similar magnitude to the upper triangle but with opposite sign.The gray boxes are associated with ecosystem transitions that we do not have reliable estimates of the slopes.This figure gives more insight into specific types of LCC on these biophysical variables.

Figure 5 .
Figure 5.The estimated sensitivities and associated uncertainty of LST to significant ecosystem transitions for (A) winter, (B) spring, (C) summer, and (D) fall seasons.The gray boxes are associated with ecosystem transitions that we do not have reliable estimates of the slopes.An interesting observation is the contrasting response of LST due to some LCC between the summer and winter seasons.
Thus, the next step is to calculate ∆X NV .Here we explain our step-by-step calculation of ∆X NV .We start with a given 1000 × 1000 m biophysical variable for 2003 and 2013.Additionally, we start with 30 × 30 m land cover maps for the years 2003 and 2013, and thus each 1000 × 1000 m pixel of a given biophysical variable includes ∼950 pixels of land cover data.We further use information from neighboring (50 km radius) pixels of the given biophysical variable based on a 50 km radius (i.e.50-pixel radius).We use these input data to calculate ∆X NV and ∆X LCC following the below steps: can calculate an LCC table from 2003 to 2013.The rows on this table are land cover in 2003 and columns are land cover in 2013.(b) From this table, we extract the percent cover change of LC1 to LC2.This value with associated ∆X LCC then becomes a data point in our regression analyses.We then have n number of ∆X LCC , hence we have n tables and n data points in our regression analysis.(c) Next, we aggregate the data by creating LCC bins in 0.1% intervals and we take the mean of all ∆X LCC observations in each bin.The main reason for binning data is to reduce the computation time and also reduce noise in the data.(d) Finally, we fit 10 000 weighted least square The steps are as follow:(a) Let us consider a 1000 m ∆X LCC pixel that is calculated in the previous section.This pixel includes ∼950 land cover pixels.Using these land cover data we

Table 1 .
Land surface temperature (LST), albedo, and evapotranspiration (ET) of different land cover types in the years 2003 and 2013.Values for each of these variables represent the mean ± standard deviation.The #obs is the number of observations where pixels were dominated by a specific land cover (i.e.land cover fraction >0.98).
Table 1 shows the mean values of biophysical variables for pixels dominated by different land cover types (i.e.land cover fraction >0.98).The retrieved mean values between 2003 and 2013 are very close.The highest annual LST and ET belong to DFs and the barren has the lowest LST and ET.The barren land cover also has

Table 2 .
Changes in the biophysical variables land surface temperature (LST), albedo, and evapotranspiration (ET) due to land cover change across the full Arctic Boreal Region (ABR).Values reported for each variable represent the maximum (Max), minimum (Min), mean (Mean), and standard deviation (StD).
the highest albedo.Forest (EF and DF) and fen have the lowest values of albedo.