Environmental conditions for alternative tree-cover states in high latitudes

Previous analysis of the vegetation cover from remote sensing revealed the existence of three alternative modes in the frequency distribution of boreal tree cover: a sparsely vegetated treeless state, an open woodland state, and a forest state. Identifying which are the regions subject to multimodality, and assessing which are the main factors underlying their existence, is important to project future change of natural vegetation cover and its effect on climate. We study the link between the tree-cover fraction distribution and eight globally observed environmental factors: mean annual rainfall, mean minimum temperature, growing degree days above 0 C, permafrost distribution, mean spring soil moisture, wildfire occurrence frequency, soil texture, and mean thawing depth. Through the use of generalised additive models, conditional histograms, and phase-space analysis, we find that environmental conditions exert a strong control over the tree-cover distribution, uniquely determining its state among the three dominant modes in ∼ 95 % of the cases. Additionally, we find that the link between individual environmental variables and tree cover is different within the four boreal regions considered here, namely eastern North Eurasia, western North Eurasia, eastern North America, and western North America. Furthermore, using a classification based on rainfall, minimum temperatures, permafrost distribution, soil moisture, wildfire frequency, and soil texture, we show the location of areas with potentially alternative treecover states under the same environmental conditions in the boreal region. These areas, although encompassing a minor fraction of the boreal area (∼ 5 %), correspond to possible transition zones with a reduced resilience to disturbances. Hence, they are of interest for a more detailed analysis of land–atmosphere interactions.

. Variables and datasets summary.

Data Analysis
After filtering and dividing the dataset, we confirm the multimodality of the tree cover distribution in high latitudes, as found by Scheffer et al. (2012) and in line with results from Xu et al. (2015), by optimising the fitting of different sums of Gaussian functions over the TCF distribution (not shown). Next, we group all data gridcells according to the modal peaks into three states: "treeless", where TCF is smaller than 20%, "open woodland", with TCF between 20% and 45%, and "forest", where 5 TCF is greater than 45%. The ensuing data analysis is aimed at two main purposes: to ascertain the impact of EVs on the TCF, and to assess whether different vegetation states can be found under the same set of EVs.
First, we evaluate the impact of the eight environmental factors on the TCF distribution using Generalised Additive Models (Miller et al., 2007). GAMs are data-driven statistical models able to handle non-linear data structures (Hastie andTibshirani, 1986, 1990;Clark, 2013); their purpose is to ascertain the contributions and roles of the different variables, thus allowing 10 a better understanding of the systems (Guisan et al., 2002). Each GAM test provides an estimate of the proportion of TCF distribution that can be explained through a smooth of one or more EVs (Staver et al., 2011b) -for instance, the formula T CF = s 1 (M T min) + s 2 (M AR) is used to assess the contribution of minimum temperature and precipitation on the tree cover fraction distribution. For each region, we repeatedly apply GAMs including different combinations of variables, and -to determine whether the sample size influences the results -we use in turn either multiple random samples of 500 gridcells each, 15 multiple random samples of 1000 gridcells each, or all the gridcells.
Subsequently, we analyse the conditional 2-dimensional phase-space between the EVs to visualise whether intersections of vegetation states in each phase-space are possible or not. To do so, we perform a kernel density estimation (KDE) of the joint distribution between the two EVs, conditioned to whether ot not the corresponding data belong to the treeless, open woodland, or forest state, and we plot the KDE together with the EVs histograms. Kernel density estimates are used to approximate the 20 probability density function underlying a set of data (Silverman, 1981(Silverman, , 1986).
Next we look at the 6-dimensional phase-space formed by mean annual rainfall, mean spring soil moisture, mean minimum temperature, permafrost distribution, wildfire frequency, and soil texture, and we divide it into classes in the following manner.
First, for every region, we divide the domain of each EV into bins. To do so, we compute the 10th and 90th percentile of the three vegetation states with respect to MTmin, MSSM, MAR, PZI, and FF. Then, for the same EVs, we select the second 25 lowest 10th and second highest 90th percentiles; these two values are the boundaries of the first and last bin, while the range in between them is equally divided into bins: 5 for MTmin, MSSM, and MAR, and 3 for FF and PZI, as exemplified in Fig. 1   10th percentile and second highest 90th percentile of the three vegetation states, with respect to the EV in use, having in mind that only one vegetation state is generally found below or above this thresholds, respectively. The remaining space is subdivided uniformly.
Afterwards, to assess our research question, we associate every gridcell of the boreal area with its vegetation state and with the class corresponding to its EVs values. Subsequently, we select two types of areas of interest, that correspond to possible alternative states: equivalent tree cover states, defined as gridcells with different vegetation state but same EVs class, e.g., an open woodland gridcell and a forest gridcell, where all the environmental variables are in the same bins;

5
fire disturbed (FD) tree cover states, defined as gridcells with different vegetation state, where the EVs are in the same bins, except for wildfire frequency, e.g., a forest gridcell with low fire frequency and an open woodland gridcell with higher fire frequency but with the remaining EVs in the same bins.
Within this last step, to take into account internal variability and the continuous evolution of the ecosystem, we consider only environmental classes that appear significantly, with a frequency above certain thresholds (generally at least 1% of the gridcells 10 with the same vegetation state). Furthermore, we test the multimodality of the TCF distribution over gridcells with equivalent and fire disturbed tree cover states. To assess this, we employ the Silverman's test against the hypothesis of unimodality (Silverman, 1981;Hall and York, 2001). Finally, to ascertain the variability of the ecosystem, we compute the standard deviation of the TCF distribution for the period 2001-2010 over the same alternative states gridcells, and we compare it with the distributions of the alternative states.

15
The entire analysis is carried out using Python 2.7.10, IPython 4.0.1, and RStudio 0.99.441.

GAM Results
A summary of GAMs results using random samples of 1000 gridcells is reported in Table 2. GAMs analysis using all the gridcells or random samples of 1000 gridcells yields similar results, with explained deviances for the former case in between the extremes of the latter, and always with statistical p-value< 0.0001. On the other hand, using samples of 500 gridcells 5 can increase the explained portion of TCF distribution at the expenses of statistical significance, due to higher p-values, and larger-scale applicability. The impact of EVs on the TCF distribution depends on the region of interest. Furthermore, the percentage of explained TCF distribution is reduced (∼ 40% maximum combined deviance explained) if we perform the analysis on broader regions than the ones here considered, i.e., on the entire boreal area at once or on the single continents. We hypothesise this is caused by 10 the different species distribution across the regions. For instance, North America is mainly dominated by "fire embracer trees", promoting high-intensity crown fires, whereas Eurasia is populated by "resister trees" in its driest regions, i.e., Eastern North Eurasia, where only surface fires are common, and fire avoiders in Western North Eurasia, which burn less frequently due to the wetter climate of this region (Wirth, 2005;Rogers et al., 2015). As a result, despite the environmental variables having We hypothesise the unexplained percentage of TCF distribution can be linked mainly to three possible causes. First, missing factors in the evaluation, for instance insect outbreaks, which are linked to climate and play an important role in the boreal forest dynamic (Bonan and Shugart, 1989), or grazing from animals (Wal, 2006;Olofsson et al., 2010). Second, deficiencies 20 in the datasets used, such as the underestimation of fire events in the boreal region (Mangeon et al., 2015), and the limited timespan of satellite observations, as fire return intervals in high latitudes can exceed 200 years (Wirth, 2005). Third, as shown later, the presence of areas where the system is in different alternative stable states under the same environmental conditions.

Phase-space results
Phase space plots for MAR versus MTmin, and MSSM versus GDD0 in Eastern North Eurasia are shown in Fig. 2. Combining we can isolate fire disturbed tree cover states, where wildfires played a major role in the timespan covered by the satellite 10 observations. A summary of the possible vegetation states found in the system is provided in Table 3. Equivalent tree cover states gridcells and fire disturbed tree cover states gridcells are represented in Fig. 3  in Table 4. Qualitative indexes for the EVs, except for ST and PZI, represent the bin into which the variable's value falls in the classification, as described in Section 2.2; the order is: very low, low, medium-low, medium-high, high, very high. Soil texture is described as belonging to the sand, loam, or clay group. Permafrost is described as sparse, discontinuous, frequent, or continuous. Each EV class contains two possible vegetation states, e.g., forest and open woodland, that are consistently found under the same specified environmental regimes.   Table 4. Classes related to equivalent tree cover states and fire disturbed (FD) tree cover states. The qualitative marks are always relative to the extremes of the EVs distributions in the region of interest, and represent the bins into which the phase-space is subdivided. percentage points (note that TCF is measured as a percentage), whereas the average standard deviation for the alternative states gridcells is 5.77 percentage points, with only one gridcell possessing a variability greater than 18 percentage points, as shown in Fig. 5. Henceforth, the bimodality of the alternative states distributions is not influenced by the variability of the TCF.  Fire is an intrinsic component of the boreal ecosystem, with an essential role in biodiversity maintenance and stand suc-10 cession. At the same time, wildfires can contribute to alternative tree cover states as a disturbing agent. In regions of reduced resilience, in fact, fires can shift the tree cover from one vegetation state to another under the same EVs. Hence, we hypothesise that a strong disturbance could permanently change the state of the ecosystem, by the combined effect of a shift in tree cover and its potential feedbacks on the environment.
We find that regions with possible alternative tree cover states encompass only a small percentage of the boreal area 15 (∼5%). However, since temperature and temperature-related environmental variables -such as permafrost distribution -exert the strongest control on the tree cover distribution and its modes, temperature changes can greatly affect forest resilience.
Therefore, under a changing climate, regions allowing for alternative tree cover states could expand both in frequency and extent.
In the context of climate change, a gradual expansion of transition zones with reduced resilience could lead to regional 20 ecosystems shifts. As climate in the boreal area is related to tree cover through numerous biogeophysical feedbacks, such as changes in albedo and transpiration, these shifts could have a significant impact not only on the structure and functioning of the boreal forest, but also on its climate.