Improving the representation of high-latitude vegetation distribution in dynamic global vegetation models

Vegetation is an important component in global ecosystems, affecting the physical, hydrological and biogeochemical properties of the land surface. Accordingly, the way vegetation is parameterized strongly influences predictions of future climate by Earth system models. To capture future spatial and temporal changes in vegetation cover and its feedbacks to the climate system, dynamic global vegetation models (DGVMs) are included as important components of land surface models. Variation in the predicted vegetation cover from DGVMs therefore has large impacts on modelled radiative and non-radiative properties, especially over highlatitude regions. DGVMs are mostly evaluated by remotely sensed products and less often by other vegetation products or by in situ field observations. In this study, we evaluate the performance of three methods for spatial representation of present-day vegetation cover with respect to prediction of plant functional type (PFT) profiles – one based upon distribution models (DMs), one that uses a remote sensing (RS) dataset and a DGVM (CLM4.5BGCDV; Community Land Model 4.5 Bio-Geo-Chemical cycles and Dynamical Vegetation). While DGVMs predict PFT profiles based on physiological and ecological processes, a DM relies on statistical correlations between a set of predictors and the modelled target, and the RS dataset is based on classification of spectral reflectance patterns of satellite images. PFT profiles obtained from an independently collected field-based vegetation dataset from Norway were used for the evaluation. We found that RS-based PFT profiles matched the reference dataset best, closely followed by DM, whereas predictions from DGVMs often deviated strongly from the reference. DGVM predictions overestimated the area covered by boreal needleleaf evergreen trees and bare ground at the expense of boreal broadleaf deciduous trees and shrubs. Based on environmental predictors identified by DM as important, three new environmental variables (e.g. minimum temperature in May, snow water equivalent in October and precipitation seasonality) were selected as the threshold for the establishment of these high-latitude PFTs. We performed a series of sensitivity experiments to investigate if these thresholds improve the performance of the DGVM method. Based on our results, we suggest implementation of one of these novel PFT-specific thresholds (i.e. precipitation seasonality) in the DGVM method. The results highlight the potential of using PFT-specific thresholds obtained by DM in development of DGVMs in broader regions. Also, we emphasize the potential of establishing DMs as a reliable method for providing PFT distributions for evaluation of DGVMs alongside RS. Published by Copernicus Publications on behalf of the European Geosciences Union. 96 P. Horvath et al.: Improving the representation of high-latitude vegetation distribution


Abstract. Vegetation is an important component in global ecosystems, affecting the physical, hydrological and 19
biogeochemical properties of the land surface. Accordingly, the way vegetation is parameterised strongly 20 influences predictions of future climate by Earth system models. To capture future spatial and temporal changes 21 in vegetation cover and its feedbacks to the climate system, dynamic global vegetation models (DGVM) are 22 included as important components of land surface models. Variation in the predicted vegetation cover from 23 DGVMs therefore has large impacts on modelled radiative and non-radiative properties, especially over high-24 latitude regions. DGVMs are mostly evaluated by remotely sensed products, less often by other vegetation 25 products or by in-situ field observations. In this study, we evaluate the performance of three methods for spatial 26 representation of present-day vegetation cover with respect to prediction of plant functional type (PFT) profiles -27 one based upon distribution models (DM), one that uses a remote sensing (RS) dataset and a DGVM 28 (CLM4.5BGCDV). While DGVMs predicts PFT profiles based on physiological and ecological processes, DM 29 relies on statistical correlations between a set of predictors and the modelled target, and the RS dataset is based on 30 classification of spectral reflectance patterns of satellite images. PFT profiles obtained from an independently 31 collected field-based vegetation dataset from Norway were used for the evaluation. We found that RS-based PFT 32 profiles matched the reference dataset best, closely followed by DM, whereas predictions from DGVM often 33 deviated strongly from the reference. DGVM predictions overestimated the area covered by boreal needleleaf 34 evergreen trees and bare ground at the expense of boreal broadleaf deciduous trees and shrubs. Based on 35 environmental predictors identified by DM as important, three new environmental variables (e.g. minimum 36 temperature in May, snow water equivalent in October and precipitation seasonality) were selected as the threshold 37 for the establishment of these high-latitude PFTs. We performed a series of sensitivity experiments to investigate 38 if these thresholds improve the performance of the DGVM. Based on our results, we suggest implementation of 39 one of these novel PFT-specific thresholds (i.e., precipitation seasonality) in the DGVM. The results highlight the 40 potential of using PFT-specific thresholds obtained by DM in development of DGVMs in broader regions. Also, 41 we emphasize the potential of establishing DM as a reliable method for providing PFT distributions for evaluation 42 of DGVMs alongside RS. Among the less explored methods to generate wall-to-wall vegetation cover predictions is distribution modelling. 82 Distribution models (DMs) are most often used to predict the distribution of a target, by establishment of statistical 83 relationship between the target (response) and the environment (predictors) (e.g. Halvorsen, 2012). The most 84 common use of DM in ecology is for prediction of species distributions (Henderson et al., 2014), but DM methods 85 have proved valuable also for prediction of targets at higher levels of bio-, geo-or eco-diversity (i.e. vegetation 86 types and land-cover types) ( The study area covers mainland Norway, spanning latitudes from 57°57'N to 71°11'N and longitudes from 4°29'E 108 to 31°10'E. Norway is characterized by a gradient from a rugged terrain with deep valleys and fjords in the western, 109 oceanic parts to gently undulating hills and shallow valleys in the central and eastern, more continental parts. 110 Temperature and precipitation show considerable variation with latitude, distance from the coast and altitude 111 (Førland, 1979). While the mean annual precipitation ranges from 278 mm in the central inland of S Norway to 112 more than 5000 mm in mid-fjord regions along the western coast, the yearly mean temperature ranges from 7°C 113 in the southwestern lowlands to -4°C in the high mountains (Hanssen-Bauer et al., 2017). 114 The vegetation of Norway is structured along two main bioclimatic gradients (Fig. 1); one related to 115 temperature/growing-season length and one to humidity/oceanity (Bakkestuen et al., 2008). Broadleaf deciduous 116 forests, regularly found in the southern and southwestern parts (the boreonemoral bioclimatic zone), are further 117 west and north (in the southern boreal zone) restricted to locally warm sites (Moen, 1999). With declining 118 temperatures northwards and towards higher altitudes, evergreen coniferous boreal forests dominate in the 119 southern and middle boreal zones. In the northern boreal zone, the coniferous boreal forests pass gradually into 120 subalpine birch forests, which form the tree line in Norway. A total of about 38% of mainland Norway is covered 121 by forests, and about 37% of the land is situated above the forest line (of which two thirds is covered by alpine 122 mountain heaths). Wetlands cover approximately 9% and broadleaf deciduous forests about 0.4% of the land area 123 (Bryn et al., 2018). 124

2.2
The AR reference dataset 125 Data obtained by in-situ field mapping, which is considered among the most reliable sources of land-cover 126 information (Alexander and Millington, 2000), is practically and economically impossible to obtain in a wall-to-127 wall format for large land areas such as countries (Ullerud et al., 2020). As an alternative, area-frame surveys 128 based upon stratified statistical sampling may provide accurate, area-representative, homogeneous and unbiased 129 land-cover and land-use data for large areas. To evaluate the three methods for representing vegetation addressed 130 in this study, we used the 'Norwegian land cover and land resource survey of the outfields' (Arealregnskap for 131 utmark) dataset (Strand, 2013), a Norwegian implementation of the mapping program LUCAS (Eurostat, 2003).

Study plots 139
Twenty out of the 1081 rectangular AR plots were selected to make up our reference dataset, AR ( Fig. 1; center 140 coordinates in Table S1). The AR plots spanned elevations from 88 to 1670 m a.s.l., with mean annual temperatures 141 between -4.0°C and 7.1°C and mean annual precipitation between 466 and 2661 mm (Table S1). The gradients of acceptably well compared to the full set of 1081 AR plots (Fig. S2). Additionally, we tested the representativeness 147 across the range of variation for a third variable (precipitation seasonality) which was later selected for sensitivity 148 experiments (see further section 4). While low values of temperature and precipitation are slightly 149 underrepresented in the 20 plots, the total range of variation was well covered. None of the tests for temperature, 150 precipitation and the additional variable (precipitation seasonality) indicate that the sample of the 20 plots deviates 151 from the full set of 1081 plots. The representativeness of the 20 plots was also tested against the full dataset of 152 1081 AR plots with regard to PFT coverage (Supplement S3, Table S3), using a Chi-square test. This test showed 153 that the two datasets are not more dissimilar than expected by chance. 154  Table S1.

Methods for representing vegetation 159
In this study, we use 'plot' as a collective term for two partly overlapping spatial units: (i) the 0.9-km 2 rectangles 160 of the AR reference dataset; and (ii) the 1-km 2 quadrats with the same centerpoint as, and edges parallel to those 161 of, the AR rectangles. The latter were used for the three methods of DGVM, RS and DM (Fig. S2). 162 Representations of the present-day vegetation for each of these 20 plots were obtained by three different methods: 163 (i) as the result of single-cell DGVM simulations for each plot; (ii) inferred from an RS vegetation map of the 164 study area; and (iii) from vegetation-type DM models (Table 1). In order to make the three methods comparable, 165 vegetation was represented by plant functional type profiles (PFT profiles), obtained by a conversion scheme 166 (Table S5 and Sect. 2.5). We define a PFT profile as a thematic representation of the land surface in a given plot 167 or a group of plots, described as a vector of relative PFT abundances, i.e. values that sum up to 1. 168 The DGVM employed in this study was the CLM4.5BGCDV (hereafter referred to as DGVM), an option provided 173 in NCAR's Community Land Model version 4.5 (CLM4.5) with vegetation dynamics, plant-soil carbon/nitrogen 174 cycle, and multi-layer vertical soil enabled (Oleson et al., 2013). In DGVM, plant photosynthesis, stomatal 175 conductance, carbon/nitrogen allocation, plant phenology and multi-layer soil biogeochemistry are described in 176 accordance with default CLM4.5, while vegetation dynamics (establishment, survival, mortality and light 177 competition) are handled separately based upon simple assumptions of environmental thresholds for establishment, 178 survival and mortality of each PFT (see supplement S6) (Oleson et al., 2013). We used DGVM in the form of 179 single-cell simulations for the 20 plots with grid-cell size set to 1×1 km (Table 1)  An inspection of the choice of atmospheric forcing, by which the CORDEX data were compared with the SeNorge 186 data used for DM, showed minimal differences (Fig. S5). Only results obtained using CORDEX data are therefore 187 shown in this paper. The 30-year CORDEX data was cycled during the spin-up. A 30-year period is consistent 188 with WMO climatological normal based on the rationale that a 30 year-period is short enough to avoid large long-189 term trends while long enough to include the range of variability. Thus, the data are not de-trended or averaged. 190 The model was run with default PFT parameters (Table S7) (Table S5)

The RS method 211
As RS product we used SatVeg (Johansen, 2009), a vegetation map for Norway with 25 land-cover classes and a 212 spatial resolution (grid cell size) of 30 m (Table 1). SatVeg is obtained by a combination of unsupervised and 213 supervised classification methods, applied to Landsat 5/TM and Landsat 7/ETM+ images within the near-infrared 214 and mid-infrared spectrum covering the period 1999-2006. While with the supervised classification, training data 215 is based on well-labelled data from the study area, during the unsupervised classification the algorithm is only 216 supplied with the number of output classes without further interference of the user. Only grid cells that were within 217 each 1-km 2 plot with a majority of their area were taken into consideration for further calculations. 218

2.4.3
The DM method 219 The distribution models (DMs) for 31 vegetation types (VT) obtained by Horvath et al. (2019) using generalized 220 linear models (GLMs, with logit link and binomial errors, i.e. logistic regression), were used for this study. The 221 VT data were collected during years 2004-2014. The DMs were obtained by using wall-to-wall data for 116 222 environmental predictors from six groups (topographic, geological, proximity, climatic, snow and land cover), 223 gridded to a spatial resolution of 100×100 m (Table 1)  for details). AUC can be interpreted as the probability that the model predicts a higher suitability value for a 229 random presence grid cell than for a random absence grid cell (Fielding and Bell, 1997). A seamless vegetation 230 map (i.e. with one predicted VT for each grid cell with no overlap and no gaps) was obtained from the stack of 31 231 probability surfaces by assigning to each grid cell the VT with the highest predicted probability of occurrence 232 within that cell (Ferrier et al., 2002). Grid cells with the majority of their area within a 1-km 2 plot were used for 233 further calculations (Fig. S6). 234

Conversion to PFT profiles 235
Harmonisation of the various vegetation classification systems was accomplished by a conversion scheme that 236 represented each grid cell (RS and DM) or polygon (AR) in each of the 20 plots with one out of the six PFTs 237 recognised by DGVM (Table S5 and Fig. S6). The scheme was obtained by expert judgements and solicited by a 238 consensus process which involved ecologists participating in the AR18x18 survey as well as scientists working 239 with RS and DGVMs. 240 We used the conversion scheme of Table S5 to generate wall-to-wall PFT maps from the original RS, DM and AR 241 datasets (Table 1) . Table S5) were left out to minimise the effect of land use, which could otherwise have brought about 245 differences in PFT profiles among the compared methods. PFT profiles were obtained for each combination of 246 method and plot. To test for deviations in PFT coverage between the methods across the whole study area, 247 aggregated PFT profiles were obtained by averaging the 20 PFT profiles obtained for each method. 248

Comparison of PFT profiles 249
To examine the overall pattern across the study area and to assess the models' ability to produce overall predictions 250 of PFTs that accord with the PFTs' overall frequency (as given by the reference) aggregated PFT profiles obtained 251 by each of the DGVM, RS and DM methods were compared with the aggregated PFT profile of the AR reference 252 dataset by a chi-square test (Zuur et al., 2007). To identify strongly deviating modelling results at a plot scale, the 253 dissimilarity between PFTs profiles obtained by each of the DGVM, RS and DM methods and the PFT profile of 254 the AR dataset for each plot was calculated by using proportional dissimilarity (Czekanowski, 1909): 255 PFT profile sums to one, the index reduces to 261 The proportional dissimilarity index is appropriate for incidence data like PFT abundances, i.e. variables that take 263 zero or positive values. The index reaches a maximum value of 1 when two objects have no common presences 264 (here, PFTs present in both compared objects) and ignore joint absences (zeros). To assess the degree to which the 265 models produce pairwise similar differences, we compared the pairwise differences between the proportional 266 dissimilarity values among methods, using a Wilcoxon-Mann-Whitney paired samples test. The aggregated PFT profiles for the RS and DM datasets did not differ significantly from those of the reference 274 AR dataset according to the chi-square test, while a significant difference was found for the DGVM profiles (Table  275 2). While the proportion of grid cells attributed to the PFT boreal NET by the RS and DM methods  In accordance with results from comparisons between aggregated PFT profiles obtained by the three methods and 291 those obtained for the reference dataset, DGVM profiles for individual plots were significantly more dissimilar to 292 the AR reference than RS and DM profiles (Fig. 2) among the three methods (Fig. 2). While no dissimilarity value for RS was above 0.50, two plots (4, 19) acted as 296 strong outliers in the distribution of DM values (cf. Fig. 2). Additionally, a comparison of proportional dissimilarity 297 between pairs of methods revealed significant differences between DGVM profiles and those obtained by

304
Visual inspection of spatial patterns of PFT profile characteristics across the 20 plots suggests that the best 305 agreement among the methods was obtained for the south-eastern part of the study area, dominated by the boreal 306 NET (Fig. 3 and Table S10). Compared to the AR reference dataset, PFT profiles obtained by DGVM were 307 strongly biased: in the north (plots 17 and 18) towards boreal NET on the cost of boreal BDT, near the west coast 308

4 Sensitivity experiments and model improvement 321
We used the results of PFT profile comparisons between DGVM and the AR reference (Fig. 3) and the results 322 obtained for the DM dataset as a starting point for exploring the possible causes of the poor performance of DGVM. 323 We first identified the three most abundant PFTs (i.e. boreal NET, boreal BDT and boreal BDS) in our set of plots 324 (Table S3). Thereafter, we identified the major VTs predicted by DM in those plots that were translated into these 325 PFTs using the conversion scheme (Table S5) (pine forest, birch forest and dwarf shrub heath, respectively; Table  326 3).  Table 3. For example, in line 337 with Fig. S12, VT 2ef and its respective PFT -boreal BDS can only establish when variable swe_10 is less than 338 380mm. 339 We explored the extent to which these additional thresholds improved the performance of DGVM on the subset 340 of six plots (i.e. 1, 2, 5, 15, 17 and 18) in which the PFT profiles are most biased compared to the AR reference 341 dataset due to the overrepresentation of the boreal NEB. In total, three sensitivity experiments were carried out by 342 a stepwise process, in each step a new threshold was added cumulatively to the previous experiment (Table 3). 343 Namely, in the first sensitivity experiment (i), we added the swe_10 threshold. In the second experiment (ii), we 344 added both swe_10 and tmin_5 as the threshold. In the last experiment (iii), we added all the three novel thresholds. 345 Only the results of the third sensitivity experiment with all the three thresholds added are reported here. Results of 346 the other two experiments are summarised in Table S13. 347   The results show that while the added thresholds for swe_10 and tmin_5 had little impact on the results (Table  352 S13), the addition of the threshold for bioclim_15 (i.e., the third sensitivity experiment) largely improved the 353 performance of the DGVM on the experimental plots explored (Fig. 4). PFT profiles simulated by this experiment 354 were more similar to those of the AR reference dataset for four out of the six plots in the experimental subset (plot 355 1, 2, 5 and 15): in plots 1 and 15, Boreal NET was correctly replaced by boreal BDS; in plots 2 and 5 boreal NET 356 was replaced by boreal BDT, BDS and temperate BDT. Addition of new threshold (bioclim_15) also reduced the 357 modelled abundance of boreal NET in plots 17 and 18, but DGVM still failed to populate these plots with another 358 PFT (Fig. 4).  Table S13.

Comparison of PFT profiles 366
The maps of PFT distributions generated by DM and RS are generally similar (Fig. S9) across most of our study 367 area. This indicates that output from DM, which is rarely used for evaluating PFT distributions from DGVMs, can 368 be used for this purpose in addition to the commonly used RS-based datasets. There are, however, some differences 369 between results obtained by the two methods near the northern Norwegian coast and in the mountain areas of 370 western Norway, which will be discussed below in more details. 371 We recognise six possible explanations for the differences in PFT profiles obtained by DGVM, RS and DM for 372 the 20 plots (see Table 5), related to the following issues: (i) the conversion scheme (ref. Table S5 which covers a range of variation in the relative abundance of graminoids and, hence, shows affinity to C3 as well 389 as to BG; and the VT '8a -Damp forest', which is usually dominated by the evergreen Scots pine and converted 390 into boreal NET, but that in some instances (e.g. after clear-cutting) is dominated by deciduous trees like Betula 391 spp. and should then be converted into boreal BDT (Bryn et al., 2018). However, a close inspection of DM shows 392 that our method reproduced similar PFT profiles as the reference dataset for all plots, except two out of 20 plots 393 (the two outliers on Fig. 2, plots 4 and 19 in Fig. 3). 394 In our case, a more complicated conversion scheme is likely to be compensated for by the sub-grid complexity 395 introduced in the process by which PFT profiles are obtained. Rather than estimating a PFT profile for the 1-km 2 396 plot directly, i.e. in one operation as in DGVM, the RS-based classes and VTs are first converted into PFTs in their 397 original resolution, and then subsequently subjected to aggregation to obtain the PFT profiles. This results in a 398 sub-grid PFT heterogeneity that could otherwise be implemented by using a more complex conversion scheme. 399

What is modelled by DGVM, RS and DM 400
The methods used in this study produce different representations of the vegetated land surface in terms of actual 401 or potential natural vegetation (Table 4). In order to model future vegetation changes and feedbacks, functional 402 type-based models like DGVM implicitly address the processes that control the distribution of vegetation ( Individual models' performance might be the reason for the two plots whose PFT profiles deviate strongly from 427 the AR reference ( Fig. 2 and Fig. 3). For plot 4, the discrepancy is due to VT "1a/1b -Moss snowbed / Sedge and 428 grass snowbed", which is represented by one of the best performing among the 31 DMs. For this VT, conversion 429 scheme bias is a more likely reason for the deviant PFT profile. For plot 19, boreal BDT is modelled because the 430 VT predicted by DM is "4a -Lichen and heather birch forest". The fact that the DM for this VT is among the 431 inferior DMs (see the ranking of individual models presented in Horvath et al. (2019)) makes this explanation 432 more likely in this case. 433

Transformation of single-DM predictions into a vegetation map 434
The performance of DM on the particular plots may also be influenced by the method chosen for transforming 435 predictions from one DM for each VT into a seamless vegetation map. Assigning to each grid cell the VT with the 436 highest predicted probability of presence in that cell, which is a commonly used method for this purpose (Ferrier 437 and Guisan, 2006), favours VTs represented by good DMs. This is brought about by good DMs having a 438 distribution of predictions that is more spread out (with larger predictions for the grid cells identified as the most 439 favourable cells) than poor DMs (Halvorsen, 2012). However, since the probability of presence for each VT was 440 predicted separately for each grid cell, the probability values for every VT vary independently of the probabilities 441 for the other VTs, throughout the study area. Thus, we regard the chance that one VT consistently outperforms 442 another VT over all the grid cells to be negligible. Alternative methods for this purpose should be tested in the 443 context of DGVM evaluation. To avoid uncertainties associated with conversion between type systems and 444 perhaps even further improve the performance of DM, we recommend exploring the option of using PFTs directly 445 as targets in DM. Direct modelling of PFTs rather than taking the detour via VT models may reduce the number 446  Table 4. 455  long-term historical climate effects (such as cooler temperature in the early 20 th century) may favour boreal BDS 495 rather than boreal NET, we consider such historical effects to have only minor impact on the already large biases 496 observed in DGVM (e.g., too much boreal NET and too few BDS). We also note that DGVM used a spatially 497 coarser CORDEX reanalysis (11x11 km) to supply high temporal resolution (6-hourly) atmospheric forcing data, 498 while the climate predictors used in DM was derived from observation-based SeNorge2 dataset with 1x1 km spatial 499 resolution and daily temporal resolution. The larger biases in CORDEX reanalysis data may also contribute to the 500 large mismatch between DGVM and the reference dataset. We have compared the average annual temperature and 501 annual precipitation of the two input datasets used in DGVM and DM to look for differences (see Fig. S4). It 502 appears that precipitation estimates by CORDEX for the 20 plots were slightly higher than SeNorge estimates, the 503 converse (but less strongly) was true for temperature. The consequences of these differences in the input data 504 might be investigated in follow-up studies. 505 Despite the shortcomings discussed above, DGVM performs reasonably well for some PFTs. One example is the 506 temperate BDT, which is correctly predicted by the model to be restricted to the southern coastal plots (Bohn et 507 al., 2000;Moen, 1999). This finding suggests that some climatically driven PFTs (i.e. temperate BDT) are well 508 implemented by the existing parameters in the DGVM used in this study.  Table S5) and by other approaches to systematisation of ecodiversity (e.g. Dinerstein et al., 2017;513 Keith et al., 2020). In particular, the number of high-latitude specific PFTs is insufficient to realistically represent 514 the biodiversity of these ecoregions, as pointed out by Bjordal (2018) and Vowles & Björk (2017). Comparisons 515 between PFT profiles obtained by DGVM and profiles obtained by DM suggest specific vegetation types that need 516 to be better represented in DGVMs, either by improving an existing PFT or by adding a new PFT (e.g. dwarf 517 shrubs vs. tall shrubs; moss dominated snow-beds, wetlands, lichens). In our study, the PFT profile of DGVM is 518 represented by the six boreal PFTs, whereas the original data for RS, DM and AR include an average of 17% (ref. 519 Table S3) of the total area that are not represented by these six PFTs (classes for "Excluded" PFT category ref. 520 Table S5). This points to the missing PFTs in the classification scheme of the DGVM, but also to the challenge 521 that certain ecosystems in our study area do not have a representation in the PFT schemes of DGVM. This is 522 exemplified by wetlands; important ecosystems that are still not represented in many of the current DGVMs. This 523 is not only problematic from the perspective of land surface energy balance (Wullschleger et al., 2014), but has 524 also implications for modelling of carbon storage and cycling, and other interactions between the land surface and 525 the atmosphere (Bjordal, 2018). In line with these studies, our results demonstrate a great potential for increasing the thematic resolution of 530 DGVMs in general and not limited to the DGVM tested here in terms of developing and parameterizing new 531 specific PFTs to be representative of the high-latitude and high-altitude habitats, and also deriving parameters from 532 observations, DMs or RS products (Bjordal, 2018;Wullschleger et al., 2014), specific for the high latitudes (Druel 533 et al., 2017). 534

Sensitivity experiments 535
Adjusting DGVM parameters so that they correspond better with environmental drivers known to be functional in 536 the high-latitude PFTs has been suggested as a measure to improve the performance of DGVM (Wullschleger et 537 al., 2014). Our sensitivity experiments demonstrate that DM results can inform DGVM parameterisation based 538 upon suitability ranges of the environmental predictors recognized by DM in determining the distribution of a 539 PFT. Most notably, we recognize that the implementation of precipitation seasonality (bioclim_15 < 50) as a 540 threshold for the establishment of NET, which has not yet been used in the DGVM, improves the distribution of 541 high-latitude PFTs simulated by the DGVM. This adds to the environmental thresholds for establishment of a PFT 542 previously used in DGVMs to restrict the predicted distribution of PFTs to realistic geographic regions (Miller and 543 Smith, 2012). Even though our sensitivity experiments focus on a limited number of additional thresholds across 544 three PFTs, this approach shows promising results and is worth exploring more extensively in future studies. 545 The importance of precipitation seasonality (i.e. bioclim_15) as a critical limiting factor for the establishment of 546 boreal NET indicates that the increased seasonality impedes growth of boreal NET. While some studies have 547 emphasized the importance of seasonal distribution of rainfall on vegetation in the semi-arid areas (Zhang et al.,548 2018), the importance of this factor for high-altitude areas is less well studied (Oksanen, 1995;Sevanto et al., 549 2006). Better representation of the processes related to the response of boreal NET to water availability, especially 550 spring-drought in DGVM, also warrants further investigation. From our results for plots 17 and 18, we notice that 551 adjusting the climatic thresholds for the establishment of boreal NET does not necessarily lead to other PFTs grow. 552 Boreal BDT and BDS can establish at both plots, but their growth rates are too slow to make them occupy a large 553 area at these plots. This implies that other environmental conditions, e.g., nitrogen availability, might play a more 554 important role in limiting the growth of BDT and BDS in the tested DGVM. The biases of the DGVM in simulating 555 BDT and BDS has been widely noticed in previous studies (Castillo et al., 2012), and remains a challenge requiring 556 more investigation in the future. 557 While going into further details of which additional PFTs should be included in DGVMs and how these and other 558 PFTs should be parameterised is beyond the scope of the present paper, we emphasize the potential of using DM 559 for improving the parameters of DGVMs. More specifically, we propose more intensive exploration of DM as a 560 tool for identification of potential environmental drivers for the high-latitude PFTs, which may enhance the 561 performance of DGVMs in high-latitude ecoregions. The specific focus of our study is the boreal region, both 562 because of the importance of these ecosystems in the climate system and because of the data availability of 563 vegetation-type DM and the field-based reference dataset (AR). However, we believe that the improved DGVM 564 parameters resulting from our sensitivity experiments may be applicable to other DGVMs such as TEM and LPJ-565 GUESS (Euskirchen et al., 2009;Miller and Smith, 2012). Also, the results from this study are likely to be 566 transferable to other high-latitude areas in the circumboreal region. 567

Conclusions 568
This study demonstrates the potential of using distribution models (DM) for representing present-day vegetation 569 in evaluations of plant functional type (PFT) distributions simulated by dynamic global vegetation models 570 (DGVMs) and for improvement of specific PFT parameters within DGVMs. By identification of the main 571 differences among PFT profiles obtained by three methods (DGVM, RS and DM) in selected high-latitude plots 572 distributed across climatic gradients in Norway, we show that PFT profiles derived from DM and RS are in the 573 same range of reliability, judged by resemblance to a reference dataset (AR). Hence, we suggest that DM results 574 can be used as a complementary evaluation dataset to benchmark the present-day DGVMs. This approach is 575 recommended when high-quality RS products are not available in desired thematic resolution or when they are not 576 able to supply proxies of other properties (such as deriving parameter improvements or PFT-specific traits). 577 Comparing the twenty PFT profiles obtained by DGVM with those obtained by AR shows a large overestimation 578 by DGVM of boreal needleleaf evergreen trees (boreal NET) and bare ground at the expense of boreal broadleaf 579 deciduous trees and shrubs. This is attributed to missing processes and PFT parameterizations of high-latitude 580 PFTs in DGVM. We use DM results to identify a new PFT-specific environmental parameterprecipitation 581 seasonalitywhich, in a series of sensitivity experiments, improves the distribution of boreal NET predicted by 582 DGVM. This new PFT-specific threshold for establishment decreases the bias of boreal NET in DGVM across 583 four out of six plots and as a result, the distribution of other high-latitude PFTs is also better represented. We argue 584 that this new threshold should be transferable to other DGVMs simulating high-latitude PFTs, and that our DM-585 based approach can be well applied to other ecosystems. 586 Further development of DGVM, such as refining parameters for existing boreal PFTs and increasing the thematic 587 resolution of PFTs for boreal areas, should be strongly encouraged to achieve a more realistic simulation of the 588 distribution of vegetation by DGVM, to increase the reliability of future predictions, and the reliability of predicted 589 vegetation feedbacks in the climate system. 590 7 Acknowledgements 591 NIBIO is acknowledged for providing access to the area-frame survey AR18X18 dataset. UNINET Sigma2 is 592 acknowledged for providing computing facilities. Geir-Harald Strand is acknowledged for providing scientific 593 assistance and Michal Torma for providing technical assistance. 594 8 Data availability. 595 The model scripts for running the DGVM are available in the GitHub repository https://github.com/huitang-596 earth/Horvath_etal_BG2020, while the script used to carry out the analysis of this study is available in the GitHub 597 repository https://github.com/geco-nhm/DGVM_RS_DM_Norway. High-resolution DM-based and RS-based 598 PFT maps are available for download at the Dryad Digital Repository https://doi.org/10.5061/dryad.dfn2z34xn 599 (Fig. S9). DGVM outputs are provided in the Table S10, Table S13 and Fig. S11. 600 9 Author contributions. 601 All authors have contributed to conceptualizing the research idea. PH curated the data and was responsible for the 602 distribution modelling and for compiling and analysing the data from all methods. HT carried out the modelling 603 and sensitivity tests using the DGVM (CLM4.5-BGCDV). PH together with AB, RH and HT were responsible for 604 writing, with all authors contributing to reviewing and editing the paper. FS, AB, TKB and LMT acquired funding 605 for this research. 606