Environmental conditions for alternative tree-cover states in high latitudes

Abstract

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 ${}^{\circ }$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 $\sim$ 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 tree-cover states under the same environmental conditions in the boreal region. These areas, although encompassing a minor fraction of the boreal area ($\sim$ 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.

Introduction

Forest ecosystems are a fundamental component of the Earth, as they contribute to its biophysical and biogeochemical processes , and harbour a large proportion of global biodiversity . However, changes in species composition, structure, and function are happening in several forests around the world . These changes originate from a combination of environmental changes, such as ${\mathrm{CO}}_{\mathrm{2}}$ concentration, drought, and nitrogen deposition , and local drivers, both anthropogenic and not, such as forest management, wildfires, and grazing . Environmental and climate changes, as well as extreme events, are likely to play a more prominent role in future decades , affecting the resilience of forests – i.e. the ability to absorb disturbances maintaining similar structure and functioning – and possibly pushing them towards tipping points and alternative tree-cover states , potentially inducing ecosystem shifts .

Increasing attention has been given to the response of ecosystems to past climate changes , and to ecosystems exhibiting potential alternative tree-cover states under the same environmental conditions, as key factors to a deeper understanding of forest resilience . To this end, in this paper, we investigate the relationship between environment and remotely sensed tree-cover distribution within the boreal ecozone. Through the use of generalised additive models (GAMs), conditional histograms, and phase-space analysis, we assess whether alternative stable tree-cover states are possible in the boreal forest, and under which environmental conditions, as understanding the mechanisms underpinning them is a key point to assess future ecosystem changes .

The boreal forest is an ecosystem of key importance in the Earth system, as it encompasses almost 30 % of the global forest area and comprises about $0.74$ trillion densely distributed trees . Despite a low diversity of tree species, a boreal forest's structure and composition depend on interactions between several factors, including precipitation, air temperature, available solar radiation, nutrient availability, soil moisture, soil temperature, presence of permafrost, depth of forest floor organic layer, forest fires, and insect outbreaks . The boreal ecozone is highly sensitive to changes in climate and can affect the global climate system through its numerous feedbacks, the most important ones related to albedo changes, soil moisture recycling, and the carbon cycle . In fact, vegetation at high latitudes can influence albedo through its distribution and through its snow-masking effect, leading to warmer temperatures . During winter, a snow-covered forest has a lower albedo than snow-covered low vegetation, as tall trees mask the snow on the ground . Additionally, differences between species distributions can affect albedo in summer, as dark conifers have a lower albedo than deciduous trees or shrubs . On the other hand, during the growing season, trees induce a cooling effect due to enhanced evapotranspiration with respect to low vegetation . Finally, the boreal forest acts as a carbon sink and is responsible for an estimated $\sim$ 20 % of the world's forest total sequestered carbon . The balance between these effects determines how the boreal forest influences climate, which, in turn, affects vegetation.

Despite its multiple roles in regulating climate, the dynamics of the boreal ecosystem regarding gradual changes and critical transitions are not yet understood . In this context, multimodality of the tree-cover distribution has recently been detected within the boreal biome . An analysis of the vegetation cover from remote sensing revealed the existence of three alternative modes in the frequency distribution of boreal trees : a sparsely vegetated treeless state, an open woodland “savanna”-like state, and a forest state. In particular, it has been observed that, over a broad temperature range, these three vegetation modes coexist ; on the other hand, areas with intermediate tree cover between these distinct modes are relatively rare, suggesting that they may represent unstable temporary states . Furthermore, it has been shown that multimodality of the tree cover does not ensue from multimodality of environmental conditions, suggesting that these three modes could represent alternative stable states acting as attractors , a stable state being the state an ecosystem will return to after any small perturbation . Hence, identifying which are the regions subject to multimodality, and assessing which are the main factors underlying their existence, is important both to understand boreal forest dynamics, and to project future changes of natural vegetation cover and their effect on climate.

We do acknowledge that vegetation and climatic variables are linked through a more differentiated set of interactions than just mean annual rainfall, temperature, and forest cover. Henceforth, to improve our understanding of the boreal ecosystem dynamics, we investigate the impact of eight globally observed environmental variables (EVs) on the tree-cover fraction (TCF) distribution. To do so, we make use of satellite products spanning the time period up until 2010, incorporating both spatial and temporal information in our analysis, and taking into account the past variability of the boreal ecosystem. Furthermore, we investigate whether the three observed vegetation modes could represent alternative stable tree-cover states. To this end, we make use of GAMs, conditional histograms, phase-space analysis, and statistical tests.

In a similar fashion, it has previously been hypothesised that tropical forests and savannas can be alternative stable states under the same environmental conditions. Evidence for bistability in the tropics has been inferred through fire exclusion experiments , field observations and pollen records , mathematical models , and satellite remote sensing .

One key piece of evidence is that the tree-cover distribution in the tropics is trimodal . In fact, multimodality of the frequency distribution can be caused by the existence of alternative stable states in the system . In the case of the tropics, multimodality could be an artefact of satellite data processing ; however, it has been suggested that this issue is not of major importance . The proposed mechanism responsible for the forest–savanna bistability is a positive feedback between tree cover and fire frequency. The same mechanism has also been employed to explore the potential of multiple stable states in a global dynamic vegetation model . Per contra, it has been suggested that trimodality of the tree-cover distribution is not necessarily due to wildfires, since it can be achieved through nonlinearities in vegetation dynamics and strong climate control . The picture is far from complete, as there is evidence that other environmental factors might play a fundamental role in controlling the tree-cover distribution .

Methods and materials

Environmental variables

We study the link between the tree-cover fraction distribution of eight globally observed environmental variables (EVs): mean annual rainfall (MAR), mean minimum temperature (MTmin), growing degree days above 0 ${}^{\circ }$C (GDD0), permafrost distribution (PZI), mean spring soil moisture (MSSM), wildfire occurrence frequency (FF), soil texture (ST), and mean thawing depth (MTD). These factors are chosen based on the work of , as they represent the main drivers of the boreal forest biome. Environmental variables can be broadly grouped into temperature, water availability, and disturbances factors.

Temperature factors include mean minimum temperature, growing degree days above 0 ${}^{\circ }$C, permafrost distribution, and mean thawing depth. Soil and air temperature are two major factors responsible for boreal forest structure and dynamics . To survive frost and dessication, during winter, coniferous trees enter a period of dormancy, characterised by the suspension of growth processes and a reduction of metabolic activity . Hence, tree growth and expansion is only possible during extended periods with air temperature above 0 ${}^{\circ }$C. We use growing degree days above 0 ${}^{\circ }$C, calculated from the NCEP/NCAR reanalysis 1998–2010 , as a measure of the extent of the growing season. Growing degree days above 0 ${}^{\circ }$C (${}^{\circ }$C yr${}^{-\mathrm{1}}$), in fact, measure heat accumulation as the sum of the mean daily temperatures above 0 ${}^{\circ }$C through a year. Furthermore, low soil and air temperatures have several other important consequences. Cold air temperatures are the main regulator of the distribution of permafrost, the condition of soil when its temperature remains below 0 ${}^{\circ }$C continuously for at least 2 years. Permafrost and low soil temperatures, on the other hand, impede infiltration and regulate the release of water from the seasonal melting of the active soil layer, inhibit water uptake and root elongation, restrict nutrient availability, and slow down organic matter decomposition . To include these effects, we use the mean minimum temperature at 2 $\mathrm{m}$ (${}^{\circ }$C), and the permafrost distribution (unitless). Minimum temperatures are obtained from the NCEP/NCAR reanalysis 1998–2010 . Permafrost distribution is extracted from the Global Permafrost Zonation Index Map , which shows to what degree permafrost exists only in the most favourable conditions or nearly everywhere.

Water availability factors include mean spring soil moisture, mean annual rainfall, and soil texture. In fact, soil moisture and water availability from precipitation are also reflected in the vegetation distribution within the boreal forest biome. Due to permafrost impeding drainage, seasonal snow melt and soil thawing can guarantee a constant supply of water during the growing season . However, this can also cause severe water loss and drought damage when trees are exposed to dry winds or higher temperatures while their roots are still encased in frozen soil and cannot absorb water . At the same time, high soil moisture reduces aeration and organic matter decomposition, promoting the formation of bogs, which in turn reduce tree growth and regeneration . To incorporate water importance, we make use of three variables: mean annual rainfall ($\mathrm{mm}$ yr${}^{-\mathrm{1}}$) from the CRU TS3.22 1998–2010 dataset , mean spring soil moisture ($\mathrm{mm}$) from the CPC Soil Moisture 1998–2010 dataset , and mean thaw depth ($\mathrm{mm}$ yr${}^{-\mathrm{1}}$) from the Arctic EASE-Grid Mean Thaw Depths product . Soil water content has also another important role, as nutrient availability and microbial activities related to nutrient cycling and organic matter decomposition depend on soil water drainage . For this reason, we employ soil texture (unitless), from an improved FAO soil type dataset , to describe the type of particles forming it, and to account for nutrient cycling and availability.

Disturbances to vegetation are represented by wildfire frequency. Nutrient cycling, organic matter accumulation, soil moisture, and soil temperature are also directly affected by recurring wildfires , which, in addition, change the albedo of the land surface, thus indirectly affecting boreal air temperatures . Additionally, forest fires can influence the composition and structure of forest communities, as plant species in boreal forests have developed different species-specific traits related to fire occurrences . These adaptations generally allow trees either to survive fires, or to promote the establishment of new individuals . Different strategies lead to different fire regimes, with implications for climate feedbacks . Hence, forest fires are a critical component of the boreal forest biome, and we quantify fire frequency (fires yr${}^{-\mathrm{1}}$) in our analysis using the GFED4 burned area dataset , and the Canadian National Fire Database . A summary of the variables we use and their origin is presented in Table .

Variables and datasets summary. Percentage tree-cover fraction indicates the proportion of land per grid cell covered by trees. Mean annual rainfall corresponds to the mean cumulative precipitation in millimetres over a year. Soil moisture is measured as water height equivalents in a $1.6$ $\mathrm{m}$ soil column. Minimum temperature refers to air temperature at $\mathrm{2}$ $\mathrm{m}$ height. Permafrost zonation index shows the probability of a grid cell to have permafrost existing only in the most favourable conditions or nearly everywhere. Fire frequency is the averaged number of fire events per year. Growing degree days above $\mathrm{0}$ ${}^{\circ }$C correspond to the sum of the mean daily temperatures at $\mathrm{2}$ $\mathrm{m}$ height above $\mathrm{0}$ ${}^{\circ }$C through a year, using 6 h measurements. Soil texture describes the type of particles forming soil, ranging from sand to clay depending on the particle size. Mean thaw depth corresponds to millimetres of thawing soil during non-freezing days averaged per year. Surface elevation refers to the topographic altitude per grid cell in metres. Land cover type describes the type of vegetation and the density of the cover, independent of geo-climatic zone.

Variable	Acronym	Units	Origin	Reference
Percentage tree-cover fraction	TCF	[%]	${0.05}^{\circ }$ MODIS MOD44B V1 C5 2001–2010 product
Mean annual rainfall	MAR	[$\mathrm{mm}$ yr${}^{-\mathrm{1}}$]	CRU TS3.22 Precipitation dataset 1998–2010
Mean spring soil moisture	MSSM	[$\mathrm{mm}$]	CPC Soil Moisture dataset 1998–2010
Mean minimum 2 m temperature	MTmin	[ ${}^{\circ }$C]	NCEP/NCAR Reanalysis 1998–2010
Permafrost zonation index	PZI	[ ]	Global Permafrost Zonation Index Map
Fire frequency	FF	[$\mathrm{fires}$ yr${}^{-\mathrm{1}}$]	GFED4 burned area dataset 1996–2012;
			Canadian National Fire Database 1980–2014
Growing degree days above 0 ${}^{\circ }$C	GDD0	[${}^{\circ }$C yr${}^{-\mathrm{1}}$]	NCEP Reanalysis (NMC initialised) 1998–2010
Soil texture type	ST	[ ]	improved FAO soil type dataset
Mean thaw depth	MTD	[$\mathrm{mm}$ yr${}^{-\mathrm{1}}$]	Arctic EASE-Grid Mean Thaw Depths
Surface elevation		[$\mathrm{m}$]	Global 30-Arc-Second Elevation Dataset
Land cover type		[ ]	Global Land Cover 2000 product (GLC2000)

To describe tree cover we make use of the percentage tree-cover fraction (%) from the MODIS MOD44B V1 C5 2001–2010 product . The MODIS tree-cover dataset has certain biases and limitations: it underestimates shrubs and small woody plants, as the product was calibrated against trees above 5m tall ; it never resolves 100 % of tree cover; it is not well resolved at low tree cover ; and may not be useful for differentiating over small ranges of tree cover (less than ca. 10 %; ), as the use of classification and regression trees (CARTs) to calibrate the dataset might introduce artificial discontinuities . Regarding the particular case of the northern latitudes, an evaluation of the accuracy of the MODIS tree-cover fraction product pointed out that the dataset may not be suitable for detailed mapping and monitoring of tree cover at its finest resolution (500 m per pixel), especially for tree cover below 20 %, and that there might be a systematic bias over the Scandinavian region . To overcome these limitations, we employ MODIS VCF data at a coarser resolution (0.05${}^{\circ }$, subsequently re-projected to 0.5${}^{\circ }$), we aggregate for most of the analysis tree-cover values into three bins encompassing the 0–20, 20–45, and 45–100 % ranges, and we exclude grid cells over Scandinavia from the analysis.

In our analysis, we investigate the use of an alternative dataset for temperatures, namely the CRU TS3.22 tmn product, for the years 1998–2010 . This dataset has a finer resolution and provides a more detailed picture of the ecosystem, albeit affected by a cold bias over Canada (see CRU TS 3.22 release notes, ). Nonetheless, it shows similar patterns to the NCEP/NCAR product. The two datasets are heavily linearly correlated, although the CRU tmn product shows lower temperatures, especially over eastern North Eurasia and eastern North America. Since our analysis is independent of variable shifts, results obtained using the CRU tmn product are essentially the same (see Supplement).

All datasets are re-projected using CDO (version 1.7.0) on a regular rectangular latitude–longitude grid at 0.5${}^{\circ }$ resolution, and divided into four main areas using approximately the Canadian Shield and the Ural Mountains as middle boundaries for North America and Eurasia: western North America (45–70${}^{\circ }$ N and 100–170${}^{\circ }$ W), eastern North America (45–70${}^{\circ }$ N and 30–100${}^{\circ }$ W), western North Eurasia (50–70${}^{\circ }$ N and 33–68${}^{\circ }$ E), and eastern North Eurasia (50–70${}^{\circ }$ N and 6–170${}^{\circ }$ W). This is done in order to preserve continuity of patterns for the environmental variables and to separate areas with different characteristics, e.g. due to oceanic influence. Note that most of Europe is excluded beforehand due to the high levels of human activity and to a possible bias in MODIS data . Subsequently, data are filtered to restrict the analysis on areas with minimum anthropogenic influence and where altitude does not play a significant role . Areas to exclude are identified using the Global 30-Arc-Second Elevation dataset and the Global Land Cover 2000 product; they correspond to sites that are either bare or flooded (codes: 15 and 19–21), subject to intensive human activity (codes: 16–18 and 22), or with elevation greater than 1200 m. The resulting datasets comprise 5848 grid cells for eastern North Eurasia (EA_E), 1559 for western North Eurasia (EA_W), 1775 for eastern North America (NA_E), and 3094 for western North America (NA_W).

Within this setup, we assume that the dataset products are suitable for our investigation.

Data analysis

After filtering and dividing the dataset, we confirm the multimodality of the tree-cover distribution in high latitudes, as found by and in line with results from , by optimising the fitting of different sums of Gaussian functions over the tree-cover fraction distribution (not shown). Next, we group all data grid cells according to the modal peaks into three states: “treeless”, where tree cover is smaller than 20 %, “open woodland”, with tree cover between 20 and 45 %, and “forest”, where tree cover is greater than 45 %. The ensuing data analysis is aimed at two main purposes: to ascertain the impact of environmental variables on the tree cover, and to assess whether different vegetation states can be found under the same set of environmental variables.

First, we evaluate the link between the eight environmental factors on the tree-cover fraction distribution using GAMs . GAMs are data-driven statistical models able to handle non-linear data structures ; their purpose is to ascertain the contributions and roles of the different variables, thus allowing a better understanding of the systems . Each GAM test provides an estimate of the proportion of tree-cover fraction distribution that can be explained through a smoothing of one or more environmental variables – for instance, the formula TCF $=s_{\mathrm{1}}\text{(MTmin)}+s_{\mathrm{2}}$(MAR), with Gaussian family and identity link (see Supplement for further details on the implementation), 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 grid cells each, multiple random samples of 1000 grid cells each, or all the grid cells.

Subsequently, we analyse the conditional 2-dimensional phase space between the environmental variables 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 environmental variables, conditioned to whether or not the corresponding data belong to the treeless, open woodland, or forest state, and we plot the KDE together with the environmental variable histograms. Kernel density estimates are used to approximate the probability density function underlying a set of data .

Next, after excluding growing degree days above 0 ${}^{\circ }$C and mean thaw depth (see Sect.  for details), we look at the 6-dimensional (6-D) 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 environmental variable into bins. To do so, we compute the 10th and 90th percentile of the three vegetation states with respect to every environmental variable except soil texture. Then, for the same variables, we select the second-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.  for MTmin; ST is instead divided according to the clay, sand, and loam groups. By doing so, we separate the range of an environmental variable where overlaps between the KDEs of the vegetation states are more likely to happen, from ranges where only one vegetation state is more likely to be found (respectively the central bins and the two most external ones). Second, we consider the partition of the 6-D phase space among the environmental variables generated by the so-computed bins. Each element of this partition is defined as a class, i.e. a class is a set of bins for the environmental variables. The idea behind this analysis is to split the 6-D environmental variable space into classes where environmental variables could be considered equal for all geographical grid cells. The question, then, is whether the tree cover could be different under the same environmental conditions.

Bin division of mean minimum temperature for eastern North Eurasia. The boundaries of the first and last bins are calculated using the second-lowest 10th percentile and second-highest 90th percentile of the three vegetation states, with respect to the environmental variable in use, having in mind that only one vegetation state is generally found below or above these thresholds, respectively. The remaining space is subdivided uniformly.

[Figure omitted. See PDF]

Afterwards, to assess our research question, we associate every geographical grid cell of the boreal area with its vegetation state and with the class corresponding to its environmental variable values. Subsequently, we select two types of areas of interest that correspond to possible alternative states:

• equivalent tree-cover states, defined as grid cells with different vegetation states but the same environmental variable class, e.g. an open woodland grid cell and a forest grid cell, where all the environmental variables are in the same bins;

• fire-disturbed (FD) tree-cover states, defined as grid cells with different vegetation state, where the environmental variables are in the same bins, except for wildfire frequency, e.g. a forest grid cell with low fire frequency and an open woodland grid cell with higher fire frequency but with the remaining environmental variables 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, i.e. with a number of grid cells per vegetation state greater than 1 % of the total amount of grid cells for that same vegetation state within the entire region (see Supplement for further details). Furthermore, we test whether the tree-cover fraction distribution over grid cells with equivalent and fire-disturbed tree-cover states is multimodal or unimodal. To assess this, we employ the Silverman's test against the hypothesis of unimodality . Finally, to ascertain that results cannot be explained by the internal variability of the ecosystem alone, we compute the standard deviation of the tree-cover fraction distribution for the period 2001–2010 over the same alternative-state grid cells, and we compare it with the distributions of the alternative states.

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

Results

GAM results

Eastern North America is the region with the highest GAM results, with more than 80 % of the total deviance of tree cover explained, and every variable except fire frequency yielding higher results than in the other three regions. Additionally, the impact of environmental variables on the tree-cover fraction distribution depends on the region of interest, as can be seen in Table . For instance, soil texture influence ranges from 9–15 to 42–52 % in western and eastern North America, respectively. A summary of GAM results using random samples of 1000 grid cells is reported in Table .

Summary of GAMs performed using random samples of 1000 grid cells each. The ranges represent the spread of results obtained with different samples, whereas the values in parenthesis correspond to the average from the samples. Statistical $p$ values are $<0.0001$ for every case. Percentages of explained deviance are a measure of the goodness of fit of each GAM . Reported values are related to the influence on tree-cover fraction distribution of mean annual rainfall (MAR), mean minimum temperature (MTmin), growing degree days above 0 ${}^{\circ }$C (GDD0), permafrost distribution (PZI), mean spring soil moisture (MSSM), wildfire occurrence frequency (FF), soil texture (ST), mean thawing depth (MTD). Values are divided within the four regions of interest, namely, eastern North Eurasia (EA_E), western North Eurasia (EA_W), eastern North America (NA_E), and western North America (NA_W).

	Deviance of TCF Explained – %
Variables	EA_E	EA_W	NA_E	NA_W
MAR	24–30 (27)	28–38 (32)	51–57 (55)	28–36 (32)
MSSM	12–20 (16)	20–29 (25)	43–53 (47)	11–21 (15)
MTmin	36–44 (40)	23–31 (27)	70–75 (72)	36–43 (40)
PZI	38–45 (42)	10–17 (13)	69–75 (71)	31–37 (34)
FF	2–9 (5)	15–20 (18)	8–13 (11)	11–19 (14)
GDD0	49–57 (54)	40–51 (46)	70–74 (71)	24–34 (28)
ST	9–18 (12)	26–35 (30)	42–52 (47)	9–15 (12)
MTD	21–33 (26)	27–37 (32)	39–46 (43)	18–30 (23)
MAR+MSSM	26–31 (28)	29–41 (34)	56–62 (59)	31–38 (34)
MTmin+GDD0	53–60 (56)	43–54 (49)	73–77 (75)	42–50 (46)
PZI+FF	42–48 (46)	34–42 (36)	70–76 (73)	34–42 (38)
All	60–67 (63)	52–58 (55)	80–85 (82)	59–65 (62)

Growing degree days above 0 ${}^{\circ }$C and mean minimum temperature are the environmental variables with the greatest influence on the tree-cover distribution, with a combined effect ranging from 42 to 77 %, in line with literature, as temperature is the main limiting factor for boreal forest . The next environmental variable in order of importance is permafrost distribution, with an impact ranging from 10–17 to 69–75 % depending on the southern extent of continuous permafrost. Water availability, as expressed through the combined effect of rainfall and soil moisture, explains 26 to 62 % of the tree-cover distribution. The two variables have a similar influence when considered alone, although MAR always has a greater impact. The impact of wildfires depends heavily on the region of interest, with FF contributing the lowest in eastern North Eurasia and the highest in western North Eurasia: 2–9 and 15–20 % respectively. Soil-related variables, namely soil texture and thaw depth, have a similar impact, generally around 30 %.

The environmental variables are not independent of each other, and hence the combined impact of multiple variables does not correspond to the sum of the single terms. For instance, PZI, MTmin, and GDD0 are highly correlated, and their combined effect is only slightly greater than the effect of each factor alone. Overall, the combined effect of all the environmental variables contributes to 52–67 % of the tree-cover fraction distribution, with the exception of eastern North America, where the cold temperatures, permafrost distribution, and rainfall gradients clearly dominate the tree-cover distribution and make up for almost 80 % of it (omitted from Table ). We obtain similar results when combining temperature-related environmental variables (GDD0, MTmin) with water-related ones (MAR, MSSM).

Performing GAM analysis using all the grid cells or random samples of 1000 grid cells 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 grid cells can increase the explained portion of TCF distribution at the expense of statistical significance, due to higher $p$ values, and larger-scale applicability. Furthermore, the percentage of explained tree-cover fraction distribution is reduced ($\sim 40\mathit{\%}$ maximum combined deviance explained) if we perform the analysis on broader regions than the ones considered here, i.e. on the entire boreal area at once or on the individual continents.

Phase-space analysis

Combining together environmental variables in phase space and performing a kernel density estimation of the joint distribution between the two environmental variables, conditioned to whether or not the corresponding data belong to the treeless, open woodland, or forest state, it is possible to locate peaks in the distributions of the vegetation states.

In many phase-space regions, environmental conditions support only a single “dominant” vegetation state. For instance, low values of GDD0 clearly denote a peak in the distribution of the treeless state. Unfortunately, GDD0 does generally not separate well between the vegetation states in the central area of its distribution, and even when combined with other variables, a clear picture does not emerge. For this reason, and for its high correlation with MTmin (Pearson's correlation coefficient $0.78<r<0.94$), GDD0 is not used in the classification. Similarly, MTD is also excluded. Nonetheless, peaks of the KDEs are not always completely disjoint, and it is possible to find intersections between the KDEs of the different vegetation states, as for the case of mean annual rainfall and mean minimum temperature with values around $400$ $\mathrm{mm}$ and $-\mathrm{7}$ ${}^{\circ }$C, respectively, where both forest and open woodland are possible. This means that the same environmental conditions can lead to different vegetation states, hinting at possible alternative states.

As a representative case, phase-space plots for eastern North Eurasia are shown in Fig. . Particularly, Fig. a represents the KDE of the joint distribution between MAR and MTmin. Each colour is associated with a vegetation state: green for forest, orange for open woodland, and purple for treeless. The isolines describe the probability of finding the three vegetation states under the specific environmental variables regimes, with intense colours indicating higher probabilities. The marginal distributions are reported on the sides of the plot in the form of histograms. The intersections of isolines marked in Fig. a show phase-space regions where the same environmental conditions can lead to different vegetation states. Similarly, Fig. b represents the KDE of the joint distribution between MSSM and GDD0, with highlighted areas where a single dominant vegetation state is supported by the environmental variables.

Results vary by region, and a complete description of all the combinations between variables is beyond the scope of this paper. Suffice to say that extremes in the distributions of environmental variables are generally associated with a single vegetation state, as in Fig. b, whereas intermediate values allow for both single states and intersections (Fig. a and b, respectively). However, these intersections consider only two environmental variables at a time and they provide only part of the total picture. Results from the classification described and discussed in Sects.  and  cover all the environmental variables at once.

Representation of the KDEs of the three vegetation states in the phase space generated by mean minimum temperature and mean annual rainfall (a), and mean spring soil moisture and growing degree days above 0 ${}^{\circ }$C (b), for eastern North Eurasia. Vegetation states are colour-coded as follows: green for forest, orange for open woodland, and purple for treeless. The isolines describe the probability of finding the three vegetation states under the specified environmental variables regimes, with intense colours indicating higher probabilities. Highlighted intersections in phase space represent areas with different vegetation states under the same environmental conditions (a), whereas the marked areas with only one dominant state hint at the unimodality of the underlying distribution (b). Marginal distributions for the variables are reported to the sides of the plots in the form of histograms.

[Figure omitted. See PDF]

6-D phase-space classification

Associating to every grid cell a class based on the values of the environmental variables reveals that in most cases (2527 classes out of 2546) there is a uniquely determined vegetation state for every class of environmental variables. However, 14 classes allow for different vegetation states, namely either treeless and open woodland, or forest and open woodland. Grid cells belonging to these classes are called equivalent tree-cover states. Furthermore, by selecting grid cells corresponding to classes differing only in the fire regime, we can isolate fire-disturbed tree-cover states, where wildfires played a major role in the timespan covered by the satellite observations (five classes). A summary of the possible vegetation states found in the system is provided in Table , divided into unimodal, multimodal, and fire-disturbed states. Equivalent tree-cover state grid cells and fire-disturbed tree-cover state grid cells are represented in Fig.  and they cover approximately 5 % of the total boreal area. Specifically, each class contains on average $29$ grid cells. Note that we excluded classes containing less than $\mathrm{1}\mathit{\%}$ of the grid cells corresponding to each vegetation state. Equivalent tree-cover states can be found in every region, with a total of 14 different environmental variable classes related to them, whereas fire-disturbed (FD) states appear consistently only in eastern North Eurasia, and consist of 5 environmental variable classes, of which 4 are also related to equivalent tree-cover states. All 19 classes are reported in Table . Qualitative indexes for the environmental variables, except for ST and PZI, represent the bin into which the variable's value falls in the classification, as described in Sect. ; the order is as follows: very low, low, medium-low, medium-high, high, very high. Precise values are reported in Table  (see Supplement). Soil texture is described as belonging to the sand, loam, or clay group. Permafrost is described as sparse, discontinuous, frequent, or continuous. Each environmental variable class contains two possible vegetation states, e.g. forest and open woodland, that are consistently found under the same specified environmental regimes.

Summary of possible vegetation states, divided as monostable, bistable, and fire-disturbed. Fire-disturbed states have a higher fire regime than the indicated counterpart. Treeless always refers to TCF $<$ 20 %, open woodland to 20 % $\le$ TCF $<$ 45 %, and forest to TCF $\ge$ 45 %.

Monostable	Bistable	Fire-disturbed
Treeless	Treeless – open	Open woodland –
	woodland	FD treeless
Open woodland		Forest – FD
		open woodland
Forest	Forest – open	Open Woodland –
	woodland	FD forest

Classes related to equivalent tree-cover states and fire-disturbed (FD) tree-cover states. The qualitative marks for fire frequency, mean annual rainfall, mean spring soil moisture, and mean minimum temperature are relative to the extremes of their distributions in the region of interest, and represent the bins into which the phase space is subdivided. Precise values for these bins are reported in brackets. Soil texture is described as belonging to the sand, loam, or clay group. Permafrost is described as sparse, discontinuous, frequent, or continuous. Each environmental variable class contains two possible vegetation states, e.g. forest and open woodland, that are consistently found under the same specified environmental regimes.

Region	Case and vegetation states	FF	ST	PZI	MAR	MSSM	MTmin	Grid cells
NA_W	1 Forest – Open Woodland	medium-low	loam	sparse	medium-high	medium-low	medium-high	27
		$[0.29;0.59]$			$[378;471]$	$[188;239]$	$[-3.6;-1.0]$
	2 Forest – Open Woodland	medium-low	clay	sparse	medium-high	medium-low	medium-high	44
		$[0.29;0.59]$			$[378;471]$	$[188;239]$	$[-3.6;-1.0]$
NA_E	3 Treeless – Open Woodland	very low	sand	frequent	low	high	medium-high	24
		$[\mathrm{0};0.07]$			$[535;647]$	$[427;490]$	$[-2.6;-0.65]$
	4 Treeless – Open Woodland	very low	sand	continuous	very low	medium-high	medium-low	20
		$[\mathrm{0};0.07]$			$[\mathrm{0};535]$	$[364;427]$	$[-4.6;-2.6]$
	5 Forest – Open Woodland	very low	sand	sparse	very high	very high	very high	58
		$[\mathrm{0};0.07]$			$[984;1607]$	$[490;598]$	$[1.3;5.9]$
EA_W	6 Treeless – Open Woodland	very high	loam	sparse	high	very low	very low	40
		$[0.59;3.18]$			$[615;663]$	$[99;257]$	$[-8.3;-4.8]$
	7 Forest – Open Woodland	very low	sand	sparse	medium-high	high	high	18
		$[\mathrm{0};0.26]$			$[568;615]$	$[327;361]$	$[-0.2;2.0]$
	8 Forest – Open Woodland	very low	loam	sparse	high	high	medium-low	20
		$[\mathrm{0};0.26]$			$[615;663]$	$[327;361]$	$[-4.8;-2.5]$
EA_E	9 Treeless – Open Woodland	very low	loam	frequent	medium-low	very high	high	35
		$[\mathrm{0};0.41]$			$[340;468]$	$[332;573]$	$[-4.5;-2.5]$
	10 Treeless – Open Woodland	medium-low	loam	continuous	very low	low	very low	34
		$[0.41;0.82]$			$[132;331]$	$[155;199]$	$[-17.9;-10.5]$
	11 Forest – Open Woodland	very low	loam	frequent	medium-low	medium-low	medium-low	23
		$[\mathrm{0};0.41]$			$[340;468]$	$[199;243]$	$[-8.5;-6.5]$
	12 Forest – Open Woodland	very low	loam	frequent	medium-high	very high	high	23
		$[\mathrm{0};0.41]$			$[468;537]$	$[332;573]$	$[-4.5;-2.5]$
	13 Forest – Open Woodland	medium-low	loam	frequent	medium-high	very high	high	21
		$[0.41;0.82]$			$[468;537]$	$[332;573]$	$[-4.5;-2.5]$
	14 Forest – Open Woodland	very low	loam	frequent	high	high	high	19
		$[\mathrm{0};0.41]$			$[537;606]$	$[288;332]$	$[-4.5;-2.5]$
FD EA_E	15 Open Woodland – FD Treeless	very low	loam	continuous	very low	low	very low	68
		$[\mathrm{0};0.41]$			$[132;331]$	$[155;199]$	$[-17.9;-10.5]$
	16 Open Woodland – FD Treeless	medium-low	loam	continuous	very low	low	very low	35
		$[0.41;0.82]$			$[132;331]$	$[155;199]$	$[-17.9;-10.5]$
	17 Open Woodland – FD Forest	very low	loam	frequent	medium-high	very high	high	11
		$[\mathrm{0};0.41]$			$[468;537]$	$[332;573]$	$[-4.5;-2.5]$
	18 Forest – FD Open Woodland	very low	loam	frequent	medium-low	medium-low	medium-low	11
		$[\mathrm{0};0.41]$			$[340;468]$	$[199;243]$	$[-8.5;-6.5]$
	19 Forest – FD Open Woodland	very low	loam	frequent	medium-high	very high	high	17
		$[\mathrm{0};0.41]$			$[468;537]$	$[332;573]$	$[-4.5;-2.5]$

Possible alternative tree-cover states over North America (left) and North Eurasia (right). The bottom five panels represent a close-up of the areas of interest ordered from west to east. Legend to be interpreted as follows: for every entry in the legend, the first name refers to the observed vegetation state in a specific grid cell, the second name corresponds to the possible alternative state found elsewhere under the same environmental conditions. Fire-disturbed tree-cover states are only present eastern North Eurasia.

[Figure omitted. See PDF]

Table  and Fig.  pinpoint the conditions and locations, respectively, of the possible alternative tree-cover states in the boreal area. To test whether the distributions of the possible alternative tree-cover states are multimodal, we employ the Silverman's test. Each Silverman's test assesses the hypothesis that the number of modes of the distributions of the alternative open woodland and treeless grid cells, and of the alternative open woodland and forest grid cells, is $\le$ 1. The tests show that the minimum number of modes to describe the distributions is two, for both cases, with $p$ values smaller than $0.001$ and $0.01$, respectively. Figure  shows the results of the Silverman's tests on the distributions of the possible alternative tree-cover states, confirming their bimodality, together with the respective tree-cover distributions. It is clear in Fig.  that both cases exhibit a decrease in frequency around 20 and 45 % tree cover.

Furthermore, we test whether the tree-cover modes can be a product of internal variability alone. To do so, we fit the distributions of the possible alternative tree-cover states using KDEs, we estimate the distances between the peaks of the distributions, and we compare them with the standard deviation of the tree-cover fraction distribution during the 2001–2010 time interval, as a measure of variability. The minimum distance between peaks corresponding to different vegetation states is $18.19$ percentage points (note that tree-cover fraction is measured as a percentage), whereas the average standard deviation for the alternative-state grid cells is $5.77$ percentage points, with only one grid cell possessing a variability greater than $18$ percentage points. Henceforth, the bimodality of the alternative-state distributions cannot be explained by the variability of the tree-cover fraction alone. A comparison between the distributions of the alternative tree-cover states, the estimated modal peaks, and internal variability is presented in Fig. .

Tree-cover fraction distribution over the grid cells where equivalent or fire-disturbed open woodland and treeless states are found (left), and where equivalent or fire-disturbed open woodland and forest states are found (right). For each case the Silverman's test verifies the hypothesis that the distribution is unimodal. The $p$ value is low in both cases, confirming the multimodality of the distributions.

[Figure omitted. See PDF]

The histogram shows the tree-cover fraction distributions of the possible alternative tree-cover states compared with tree-cover fraction internal variability. Purple bars refer to treeless/open woodland states, and green to forest/open woodland states. The black lines and the orange dots represent the kernel density estimate fittings of these distributions and the locations of their modal peaks, respectively. Internal variability of the tree-cover fraction distribution for the period 2001–2010, computed as the standard deviation of the distribution, is 5.77 percentage points, and is represented as the orange error-bars. The minimum distance between peaks corresponding to different vegetation states is 18.19 percentage points and is higher than what internal variability could explain.

[Figure omitted. See PDF]

Notably, equivalent tree-cover states generally fall into two categories: either they possess intermediate values for the environmental variables, or they have contrasting ones. For instance, case number 1 in Table  is characterised by medium or intermediate values for all the environmental variables, whereas case number 6 shows high values for FF and MAR, but very low for MSSM and MTmin. The first category, with intermediate values, can be associated with transition zones, when passing from an environmental variable class where only a single vegetation state is dominant, to a class where another state is dominant. As a result, the observed tree-cover fraction distribution can oscillate between the two states. The second category, on the other hand, relates to classes where at least one of the environmental variables has a value contrasting with the remaining ones. For instance, in case 8, PZI, MAR, MSSM, and ST all possess values generally associated with forest states; however, MTmin is low, preventing tree growth. This possibly creates a limit cycle where the ecosystem alternates between the different alternative states. Fire-disturbed tree-cover states, instead, can be grouped into three categories. The first category is represented by classes where the vegetation state with the lowest tree cover is disturbed by fire, and the one with highest tree cover corresponds to one of the existing equivalent tree-cover states (case 16, 18, and 19). The second category is the opposite: the vegetation state with the highest tree cover is disturbed and the one with the lowest tree cover is found among the equivalent tree-cover states (case 17). The third category corresponds to the first one, but neither of the vegetation states is found among equivalent tree-cover classes (case 15, although very similar to case 10).

Classification results suggest that environmental variables exert a strong, albeit sensitive, control over the tree-cover distribution. Depending on the conditions, only one of the three possible vegetation states is attained; for instance, in eastern North America, classes with very low mean annual rainfall and mean minimum temperature (MAR below $500$ $\mathrm{mm}$ $\mathrm{yr}$${}^{-\mathrm{1}}$ and MTmin lower than $-\mathrm{9}$ ${}^{\circ }$C; see Supplement) are associated with treeless grid cells. In $95\mathit{\%}$ of the grid cells, environmental conditions uniquely determine the vegetation state. However, in transition zones with intermediate or contrasting conditions, it is possible to find multiple vegetation states with the same environmental regimes. In such zones, disturbances could shift the system between the possible alternative states. In this sense, fire is part of the environment both as a variable and as a disturbance. Strong fire events in transition zones can determine which of two alternative states the system will fall into. On the other hand, changes due to fire events in a stable area should be reabsorbed with time, unless they are so dramatic as to produce changes in another main environmental variable, creating a new transition zone.

Discussion

The link between environmental variables and tree-cover fraction varies within the four boreal regions considered here, as described in Sect. , and is stronger in eastern North America, where the cold temperatures, permafrost distribution, and rainfall gradients dominate the tree-cover distribution. Furthermore, the percentage of explained tree-cover fraction distribution is greatly reduced when performing the analysis on broader regions, such as the entire boreal area at once or on the single continents. We hypothesise that this is caused by the different species distribution across the regions and by the different species-specific adaptations to the surrounding environment. For instance, North America is mainly dominated by “fire-embracing trees”, promoting the accumulation of fuel and the occurrence of high-intensity crown fires. On the other hand, Eurasia is populated by “fire-resistant trees” in its driest regions, i.e. eastern North Eurasia, where only surface fires are common, and “fire-avoider trees” in western North Eurasia, which burn less frequently due to the wetter climate of this region . As a result, despite the environmental variables having different distributions, the general response of the tree-cover distribution in the four regions is similar, but the impact of each individual environmental variable varies within the regions.

Minimum temperatures and growing degree days are the most influential environmental variables for the boreal tree-cover fraction distribution, as can be seen in Table . Nonetheless, their combined effect does not fully explain the tree-cover distribution, as a more diverse set of variables and feedbacks plays a role. Additionally, the environmental variables are not independent of each other, and hence the combined impact of multiple variables does not correspond to the sum of the single terms. Furthermore, the overall effect of the environmental variables is not able to fully explain the tree-cover distribution. We hypothesise this 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 , or grazing from animals . Second, deficiencies in the datasets used, such as the underestimation of fire events in the boreal region , and the limited timespan of satellite observations, as fire return intervals in high latitudes can exceed 200 years . Third, supported by the multimodality of the boreal forest and by the results presented in Sect. , the presence of areas where the system is in different alternative stable states under the same environmental conditions.

By linking tree-cover distribution to a 6-D phase space formed by environmental variables, we show that under most environmental conditions, the tree-cover fraction distribution is uniquely determined, i.e. it is unimodal, suggesting a strong control of the vegetation by means of the environment. In this sense, the three different modes of the boreal tree-cover distribution represent three distinct stable tree-cover states that do not generally appear under the same environmental conditions. However, we find areas where the tree-cover fraction distribution is bimodal under the same environmental conditions, suggesting the existence of possible alternative states, as depicted in Fig. . These areas are characterised by either intermediate or contrasting environmental conditions, possibly creating limit cycles that allow alternative tree-cover states. Furthermore, these areas seem to exhibit a reduced resilience, since disturbances, such as wildfires, appear to be able to shift the vegetation from one state to the other, as in the case of fire-disturbed tree-cover states. Particularly, eastern North Eurasia is the region with the greatest extent of possible alternative tree-cover states, and it is the only region where fire-disturbed states are found, hinting at a greater susceptibility of its forest resilience.

Environmental conditions control the tree-cover distribution in high latitudes, pushing its vegetation towards three distinct tree-cover states. This hints at the presence of feedbacks between the vegetation and the environment, which are able to stabilise the vegetation cover in three different ways. However, the environment is influenced by the forest cover state through albedo, water evapotranspiration , and nutrient recycling. Thus, changes in climate and environmental variables will trigger feedbacks from the vegetation that can either further amplify or dampen the initial changes. In particular, areas of reduced resilience where alternative tree-cover states are found, i.e. what we call transition zones, will be affected. As the classification results suggest, environmental variables drive the ecosystem towards seemingly stable states and away from intermediate unstable ones, resulting in the multimodality of the tree cover. Thus, disturbances in transition zones could cause a rapid ecosystem shift regarding tree cover. Henceforth, it is important to better understand the interplay between environmental variables and tree cover.

Additionally, there are other factors playing a role in the dynamics of the boreal forest, both at local and larger scales. For instance, the understory vegetation acts as an important driver of soil fertility, influencing plant growth and tree seedling establishment . An increased nitrogen deposition may promote accumulation of organic matter and carbon in boreal forest . At the same time, its effects on the forest floor and soil processes might decrease forest growth . Despite its importance, there is a lack of knowledge regarding the impact of understory interactions at large spatial scales, and the contribution of climate-change drivers . For these reasons we could not take it into account in our study. Another missing factor is nitrogen (N), as plant growth in the boreal forest is thought to be generally N-limited . Additionally, herbivore grazing is also influenced by N fertilisation , with the potential to affect feedbacks involving soil nutrient cycle and plant regeneration . However, globally distributed datasets for N availability and grazing pressure suitable for our analysis are not yet available. Local topography also plays a role, as the low solar elevation angle at high latitudes accentuates the effect of ground characteristics such as slope and aspect , affecting temperature and soil moisture. Finally, micro-topography, such as shelter from boulders, can increase resistance to disturbances by creating small-scale refugia , thus locally increasing the resilience of the forest.

In the context of climate change, understanding transition zones at large scales is necessary for assessing future projections of vegetation cover. Climate change is impacting the boreal area more rapidly and intensely than other regions on Earth; for instance, surface temperature has been increasing approximately twice as fast as the global average . Temperature is a key variable in this region, as it is connected with tree growth and mortality cycles, with permafrost thawing and the hydrological cycle, and with disturbances, such as wildfires and insect outbreaks . Particularly, air and surface warming can increase the frequency and extent of severe fires , and promote more favourable conditions for insect outbreaks . At the same time, climate-change influences the resilience of boreal forest stands , making them more susceptible to abrupt shifts due to disturbances. As temperature increases and permafrost thaws, it is more likely to find intermediate conditions where alternative tree-cover states are possible. For instance, a study on the southern part of the eastern North American boreal forest has shown that an increased disturbance regime, together with the superimposition of fires and defoliating insect outbreaks, can cause a shift between alternative vegetation states . Furthermore, there is strong evidence that certain types of extreme events, mostly heatwaves and precipitation extremes, are increasing under the effect of climate change . Such events could foster areas with contrasting environmental conditions, further weakening the stability of the boreal ecosystem, and increasing its susceptibility to shifts.

Conclusions

Through the analysis of generalised additive models, we find that the environment exerts a strong control over the tree-cover distribution, forcing it into distinct tree-cover states. Nonetheless, the tree-cover state is not always uniquely determined by the variables at use. Furthermore, the response of vegetation to the environment varies in the four regions considered: eastern North America, western North America, eastern North Eurasia, and western North Eurasia.

By means of a classification, we analyse the 6-D phase space formed by mean annual rainfall, mean minimum temperature, permafrost distribution, mean spring soil moisture, wildfire occurrence frequency, and soil texture. We find several environmental conditions under which alternative tree-cover states are possible, broadly falling into two categories: with contrasting environmental features, e.g. high rainfall but low temperature, or with intermediate environmental values. In regions under these environmental conditions, the tree cover exhibits a reduced resilience, as it can shift between alternative states if subject to forcing.

As fires can shift the tree cover from one vegetation state to another in regions of reduced resilience, we find support for the hypothesis that a strong fire 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.

Finally, we find that regions with possible alternative tree-cover states encompass only a small percentage of the boreal area ($\sim$ 5 %). However, since temperature and temperature-related environmental variables exert the strongest control on the tree-cover distribution and its modes, temperature changes could greatly affect forest resilience and cause an expansion of regions with alternative tree-cover states.

In the context of climate change, a gradual expansion of transition zones with reduced resilience could lead to regional ecosystem shifts with a significant impact not only on the structure and functioning of the boreal forest, but also on its climate.

Data availability

All data, scripts, and information necessary to reproduce the authors' work are archived by the Max Planck Institute for Meteorology. They are accessible at http://hdl.handle.net/11858/00-001M-0000-002C-4F6F-3 or by contacting Carola Kauhs at [email protected].

Information about the Supplement

All data necessary to reproduce the paper are provided in the Supplement.

The Supplement related to this article is available online at doi:10.5194/bg-14-511-2017-supplement.

Both authors designed the research. Beniamino Abis performed the research. Both authors contributed to the discussion and interpretation of the results. Beniamino Abis wrote the first draft of the manuscript and both authors contributed to its final draft.

The authors declare that they have no conflict of interest.

Acknowledgements

We would like to thank the Integrated Climate Data Center (ICDC, http://icdc.cen.uni-hamburg.de/) of the University of Hamburg for helping with the retrieval of datasets. We thankfully acknowledge Gitta Lasslop for helping us review the manuscript, and Fabio Cresto Aleina for his useful comments. We also thank the International Max Planck Research School on Earth System Modelling (IMPRS-ESM) for supporting this work. Finally, we thank our editor A. Ito and three anonymous referees for their constructive comments on the paper. The article processing charges for this open-access publication were covered by the Max Planck Society. Edited by: A. Ito Reviewed by: three anonymous referees