This manuscript describes a modeling study to understand the climatic and biogeochemical controls on vegetation change in sub-arctic Sweden. The authors parameterized the LPJ-GUESS dynamic vegetation model with the principle plant types found in their study area around Abisko, Sweden. Using a very high resolution local gridded climatology and downscaled GCM output, the authors then ran LPJ-GUESS in a historical and series of future climate scenarios, and in a range of sensitivity tests controlling for different processes including CO₂ fertilization and nitrogen cycling. The authors conclude that while climate has an overarching control on vegetation composition and position of the treeline, nitrogen availability exerts a very important influence on vegetation dynamics.

In general, this study is well designed and the methodology and sensitivity tests follow generally accepted protocols for dynamic vegetation modeling experiments. The manuscript is well written and easy to follow. The presentation could be improved, in particular the figures, and I do have a few scientific comments that should be addressed, but overall this manuscript should be suitable for publication in Biogeosciences after moderate revision.

Reply:   Thank you for these encouraging words and good suggestions for improvements to the manuscript. Below follow a few discussion points to your questions and point-to-point answers, as well as an updated manuscript with your suggestions incorporated or discussed.

General comments

I would like to see more discussion of the processes that were not included in the model, concentrating on the following:

1. Seed dispersal limitations to vegetation change: A great deal of modeling work has gone in to understanding the role of seed dispersal in limiting plant migration rates, particularly along vertical gradients similar to the principle one at Abisko. The pioneer in this research has been Heike Lischke with her TreeMig Model (Lischke et al., 2006), which has been applied to the Arctic (Epstein et al., 2007), and some representation of dispersal and migration has even been incorporated into a version of LPJ developed in part by your colleagues in Lund (Lehsten et al., 2019). Of course, dispersal limitations are not the whole story of what limits plant migration (Scherrer et al., 2020). Nevertheless, some further discussion on this topic and additional citations would be welcome in this manuscript, as it is effectively missing at the moment.

Reply:   Thank you for a good comment! As you state, a great deal of research has gone into the study of this topic, where the TreeMIG-model stands out as a pioneering example. Seed dispersal limitations were only briefly touched upon in our Introduction, where we state that latitudinal treelines might be more sensitive to dispersal constraints than elevational treelines (see L49-50), as has been discussed by Rees et al. (2020). LPJ-GUESS assumes the presence of seeds in all grid cells, meaning that simulated species or PFTs can establish once the climate is favourable, as defined by each individual’s predefined bioclimatic limits (Smith et al. 2014).  While seed dispersal limitations might play a role in constraining larger scale, latitudinal, migration rates, we do not believe that it will play any major role in our case, across distances of the order of tens of metres. This is evidenced by the presence of tree saplings above the current treeline in the Abisko area (Sundqvist et al., 2008) and that some of these stunted mountain birch individuals can be a few decades old (Hofgaard et al., 2009). The latter suggesting that these seedlings have been able to spread and establish above the treeline in historic times, but not develop into mature trees. There should thus not be any constraint on the availability of the birch seeds in the seedbank above the treeline (Truong et al., 2007). According to your suggestion we added a section in the discussion (see L473-487) about this topic.

“LPJ-GUESS assumes the presence of seeds in all gridcells and PFTs may establish when the 20-year running average climate is within PFT-specific bioclimatic limits for establishment. This assumption may overlook potential constraints on plant migration rates such as seed dispersal and reproduction. On larger spatial scales, it is likely that lags in range shifts would arise from these additional constraints (Rees et al., 2020; Brown et al., 2018). Models that account for dispersal limitations generally predict slower latitudinal tree migration than models driven solely by climate (Epstein et al., 2007). However, on smaller spatial scales, the same models predict competitive interactions to be more dominant in determining species migration rates (Scherrer et al., 2020). In a seed transplant study from the Swiss alps, seed viability could not been shown to decline towards the range limits of eight European broadleaved tree species (Kollas et al., 2012; Körner et al., 2016). Similarly, gene flow above the treeline could not be shown to be limited to near-treeline trees in the Abisko region (Truong et al., 2007). Furthermore, tree saplings have been reported to be common up to 100m above the present treeline (Sundqvist et al., 2008; Hofgaard et al., 2009). As environmental conditions improve, these individuals may form the new treeline. Thus, on the scales considered in this study, we do not regard dispersal limitations as a major factor in limiting range shifts of trees.“

2. Permafrost effects: While it is mentioned that the study area lies in the permafrost zone, it is not really discussed how changing permafrost dynamics; deepening of the active layer, changes in effective rooting depth, changing water table depth, etc. affects vegetation. At the ultra-high resolution used in the model simulations, it might be important to account for how ground freezing and soil vertical and horizontal movement caused by frost (or lack of it in the future) could affect survival and competition among the various plant functional types simulated. While I appreciate that a full treatment of permafrost dynamics may be beyond the scope of the present study, it would be good to have some further discussion/speculation of how this process could influence the results and conclusions of the modeling work performed here.

Reply: Permafrost aspects are important in the Arctic. Our study area is located within a zone of discontinuous permafrost. In our study domain, the presence of permafrost is limited to wetlands, e.g., the Stordalen area, or at the highest elevations (Callaghan et al., 2013). As we do not include any analysis of wetlands in our study, the only permafrost zones in our study domain are located on the mountains above treeline. Thus, permafrost and treelines in our area are far separated and permafrost effects on the treeline or shrub growth will be of minor importance for our results. Permafrost dynamics may be important for vegetation dynamics at the highest elevations in our historic simulations.

The model simulates the effect of frozen ground and permafrost (where present) on water availability, though the effects of frost horizontal and vertical movement are not taken into account. Such processes could affect survival and competition among the plant functional types, especially in the seedling stage when plants are most vulnerable to mechanical disturbance. These effects could be relevant to treeline dynamics at the high grid resolution of our study, but are not accounted for by our model. We have added a few lines about permafrost dynamics in the results (L279-281; 348-350) and discussion (L 456-463).

L279-281

“Simulated permafrost with an active layer thickness of <1.5 m was present at elevations down to 560 m a.s.l. in a few gridcells, but was always well above the treeline. More shallow permafrost (<1 m) was only present in gridcells at elevations of 940 m a.s.l. and above.”

L 346-348

“Permafrost with an active layer thickness of <1.5m disappeared completely from our study domain in all scenarios except the coldest (GFDL-ESM2M-RCP2.6) where it occurred in a few gridcells at elevations of approximately 600 m a.s.l. However, the shallow permafrost (<1m) had disappeared also in this scenario.”

L 456-463

“Permafrost was only present at the highest elevations during the historic simulation but had disappeared from the landscape at the end of the century for all except the coldest scenario (GFDL-ESM2M-RCP2.6). The simulated permafrost was however always well above the treeline and did not have a significant impact on the treeline advancement. While some aspects of ground freezing are accounted for in the model, soil vertical and horizontal movement caused by frost, and amelioration of such effects in the warmer future climate are not. Such processes could affect survival and competition among the plant functional types, especially in the seedling stage when plants are most vulnerable to mechanical disturbance. These effects could be relevant to treeline dynamics at the high grid resolution of our study but are not accounted for by our model.”

3. Slope and aspect effects: It is mentioned that the study area is hilly or even mountainous; at the resolution of the model, how were slope and aspect handled? Particularly in a high-latitude situation with low sun angle, slope and aspect must be very important in influencing the surface radiation budget, soil temperature, and snowpack dynamics. If slope and aspect were not considered in this study, some explanation of why is required, and similarly to the point above, the authors should include some discussion of the potential effects that this could have on their results. Furthermore, as the resolution of the modeling approaches micrometerological scale, it would be helpful to have some further discussion of how the lakeshore climate may be different from areas further away, e.g., with respect to wind speed and the radiative environment.

Reply:   The questions about the influence of slope and aspect will of course be of importance when modelling on this scale. Elevation (i.e. lapse rate) is the main driver of mesoscale temperature variation in our climate dataset.  Effects of mountainside aspect are included in the climate dataset by Yang et al. (2011) as local variations to the temperature data alongside the lake effect by Lake Torneträsk.  Both surface temperature and radiation are input variables to the model. Any relationships between surface temperature and radiation are not calculated within the model but are implicitly included in the dataset. We have therefore included a statement of these properties of the climate dataset in the methods section (see L 166-167). For a full description of the climate dataset we refer to Yang et al. (2011; 2012). Soil temperature and snowpack dynamics are determined by the pre-calculated surface air temperature and precipitation inputs, and soil properties, with no explicit use of aspect or slope in the model’s process descriptions. Effects of wind is not regarded in the model and we do not use wind as input to any parts of the model.

4. Linked to the point above on slope and aspect, I would have liked to see some more discussion of the spatial heterogeneity of the snowpack. Again, at the model resolution and over the spatial domain considered, I would imagine that the formation of snowdrifts and other snowpack variation is important for soil temperature and moisture, plant survival, and N cycling. Numerous studies have demonstrated that wind and slope/aspect have a strong influence on the depth and density of snow in snowdrifts and on the rate and timing of snowmelt. This spatial heterogeneity in snowpack depth and melt rate affects winter surface temperatures and therefore survival of plants at and above the treeline, and growing season soil moisture (there are many studies on this topic but one example is Walker et al., 1999). Again, while a full treatment of snowpack heterogeneity might be beyond the scope of the study, some more discussion of this important process, and how it might influence the region around Abisko specifically, is warranted.

Reply:   Our version of the model does not include any formulation of snowdrift or wind compaction. The simulated spatial heterogeneity of the snowpack in the model is thus minimal (<1%) for the winter (DJF) months. Snow trapping processes or any potential snow-shrub feedbacks (Sturm et al., 2001) are therefore not included in our simulations. These processes are undoubtedly important for soil temperatures and subsequently mineralisation of soil organic matter. The model accounts for snow insulation effects on soil temperature, but do not have any representation of frost damage, which might affect seedling survival and be ameliorated by an insulating snow cover. As the model does not account for heterogeneity over the landscape in the compaction and drifting of snow, this might lead to an overestimation of winter soil temperatures above the treeline, where high winds and low roughness tend to deplete snow cover, and an underestimation of winter soil temperatures in the forest, where the opposite is true. We have extended the discussion to include these feedbacks (see L 443-450).

“Furthermore, our model does not include any wind related processes such as wind mediated snow transport or compaction. Thus, our simulations result in a homogenous snowpack during the winter months with no differentiation in sheltering or frost damage that may result from different snow and ice properties. Sheltered locations in the landscape are known to promote survival of tree saplings (Sundqvist et al., 2008). For N cycling this may also mean that suggested snow-shrub feedbacks (Sturm et al., 2001; Sturm, 2005) are not possible to capture with the current version of our model. While overall rates of advance were captured, local variations arising from physical barriers such as steep slopes, stony patches or anthropogenic disturbances were consequently not possible to capture as these processes are not implemented in the model.”

Any explicit representation of physical barriers, e.g., stones, or processes such as snowdrift resulting from topography are also not included in our simulations. We have updated the manuscript to clarify that these processes are not included in the discussion section (see L 446-447). The analysis of the transects are intended as a way of clarifying how these barriers in the landscape influence our simulated treeline advance.

5. Given the overall importance of N cycling for the results of this study, it would be helpful to have an overview of the N module in LPJ-GUESS. In particular, I would like to understand how biological N fixation is represented and if certain PFTs (e.g., something representing Alnus spp.) can be advantaged in nitrogen poor settings because they are capable of enhanced N fixation especially with warmer temperatures.

Reply: A more detailed description of the N cycle in LPJ-GUESS have been included in the manuscript (see L 125-128). Full details and equations, which are too extensive to repeat in the present paper, are presented in the cited paper by Smith et al. (2014).

“Biological N fixation is represented by an empirical relationship between annual evapotranspiration and nitrogen fixation (Cleveland et al., 1999) and occurs differently within each patch. Additional inputs of nitrogen to the system occur through nitrogen deposition or fertilisation. Nitrogen is lost from the system through leaching, gaseous emissions from soils and wildfires. For a full description of the nitrogen cycle in LPJ-GUESS, see Smith et al. (2014).”

Specific comments

Line 148

It is mentioned that three replicate patches in each gridcell are used for LPJ-GUESS. It is worth going in to a little more detail here to justify this choice of the number of replicate patches. As I understand, each patch in LPJ-GUESS is meant to represent an area of 0.1 ha. With a 50m grid (cells of 0.25 ha), three replicates effectively makes an explicit representation of the entire gridcell, no?

Reply:   The replicate patches are intended to give an estimation of landscape-scale heterogeneity within a gridcell or stand that might arise from spatial variation in stochastic processes and histories in the model (e.g. stochastic establishment, mortality and patch destroying disturbance events). No assumptions are made about how the patches are distributed within a wider area, they are merely a statistical sample of equally possible disturbance/demographic histories across the landscape of a grid cell. We have adjusted the text to improve clarity on lines 147-148.

“No assumptions are made about how the patches are distributed within a gridcell, they are a statistical sample of equally possible disturbance/demographic histories across the landscape of a gridcell.”

Line 153

Further to my general comment above, please explain how slope and aspect are incorporated into this gridded climatology.

Reply:      We have added a clarification of how the dataset was constructed. In general, the field measurements were set up in forms of transects with temperature loggers. These transects were selected to cover variations in mesoscale climate patterns that arise from elevation (i.e., lapse rate), and local variations (aspect and lake effects). See lines 161-168.

“The field measurements were conducted in form of transects that captured mesoscale climatic variations, i.e., lapse rates. In addition, the transects were placed to capture microclimatic effects of the nearby lake Torneträsk and variations in radiation stemming from mountainside aspect. The temperature in the lower parts of the Abisko valley in the resulting dataset was influenced by the lake with milder winters and less yearly variability. At higher elevation, the temperature was more variable over the year and the local scale variations were more dependent on the different solar angles between seasons and mountainside aspect (Yang et al., 2011; Yang et al., 2012)”

Line 157 and Fig S1.1

From looking at the figure I don’t really see how temperature is “more variable” with increase in elevation. Perhaps some descriptive statistics would be more useful here

Reply:      We acknowledge that the statement about climate variability is not evident from the figure. It is however a feature of the climate dataset described closer in Yang et al. (2012). We have updated the figure with larger elevational bands for clarity and values written out in the figure. In addition, we have added a panel with the standard deviation for each elevational band and month. This panel may be interpreted as the magnitude of the local effect at each elevational band and month.  

Line 167-169

Soil edaphic controls on vegetation are an important part of treeline and subarctic biogeography; it is even mentioned on this line. So why not make any attempt to account for spatial variation in soils? Although the spatial resolution is still a bit coarse, why not use the pedometrics-based Soilgrids250 (Hengl et al., 2017) instead of simply prescribing the same soils everywhere? Could you have done some sensitivity tests to quantify the model response of vegetation distribution and treeline to different soil physical properties?

Reply:      We agree that soil factors are undoubtedly important also for soil nutrient cycling and storage. We assumed a uniform soil texture within our study domain sourced from our standard soil dataset (Batjes, 2005) as input to the model. We came to this decision after we had contacted the Swedish Geological Survey (SGU), which we judged as the best source of more detailed soil data or surveys. A more detailed survey of soil texture for this area was however non-existent. The dataset by Hengl et al. (2017)is indeed impressive and includes a few processes known to create variations in soil texture over the landscape. However, we would argue that a dataset with higher resolution is not necessarily a more reliable or accurate dataset. Furthermore, the variations in soil textures within our study domain are not larger than approximately 3% for the clay or sand fractions. Such small variations in soil texture would not generate any significant changes to the landscape heterogeneity. We did some sensitivity tests in a sub-section of our domain during the preparation phase of the simulations. While some factors such as soil organic carbon content was affected, more drastic variations in soil texture are needed to affect vegetation distribution or treeline dynamics to any large degree. We have not included these sensitivity tests in our results. We have added a sentence in the discussion (L451-454) about this limitation in soil texture heterogeneity in our model input.

 “High-resolution, local observations of vertically-resolved soil texture and soil organic matter content (see, e.g. Hengl et al. (2017) for an example compiled using machine-learning) have the potential to improve the spatial variability of modelled soil temperatures and nutrient cycling in our study domain. We will investigate this uncertainty in future studies.”

Line 418

Where are these transects? Call out the supplementary figure here or even better refer to an overview map (see comment below). How were the locations and orientations of these transects chosen? 

Reply: The figure has been referenced in the text on line 235. We also include a short sentence in the methods section (see L 235-237) and figure caption (supplementary materials) about how these transects were selected.

Comments on presentation

I would appreciate seeing an overview map or aerial photo of the study area showing topography and the location of the lake, any rivers, and settlements, roads, etc. I would also like to see at least an inset map showing the location of the study area within Europe and Sweden.

Reply:      We have provided an inset map of Sweden and the location of Abisko in figure 2. In this figure we added contour lines to mark the landscape topography. We also marked the lake and the location of the Abisko scientific station (ANS) in this figure. We thank the reviewer for these suggestions which we believe improves the maps a lot. We do not include settlements, roads and rivers as our simulations are not affected by these features.

Fig. 2

What is the white area in these maps? Why are the colors used for the PFTs in Fig. 2a not the same as those used in Fig. 3? Please harmonize. Please add a scale bar to these figures, and perhaps one or two longitude and latitude tick marks/labels. As many readers look only at the figures, or the figures first, it would be helpful to spell out the PFT names in the figure legend here instead of making the reader refer back to an additional table or text to decode these.

Reply:      We have updated the figure and figure caption in accordance with these suggestions.

Fig. 3

Harmonize the colors with Fig. 2a

Reply:   We have updated the figure with harmonized colors as in figure 2.

Fig. S1.1

The figure caption appears to be cut off

Reply:   Fixed.

Supplement S1 Table 1

What is reported in the column “Reported (van Bogart et al. 2011)”? What are these units of?

 Reply:     These are estimated treeline migration in elevational meters reported by van Bogart et al (2011). We have updated this table to clarify this.

Fig. S1.4

What is the gray scape in plotted in the background of the map? What is the white area?

 Reply:     The gray scape in this figure is the landscape relief and the white area is Lake Torneträsk. We have clarified this in the figure caption.

Fig. S1.6e

What is the principle control on annual shortwave radiation? Is it cloud cover? This could also be discussed in the main text.

Reply:      The principal control of the annual shortwave radiation within the global climate models (GCMs) is cloud cover. However, we use the monthly bias adjusted shortwave radiation-output from the GCMs as input to our model.  The annual values are averages of the monthly output provided by the GCMs. We have updated the figure caption to clarify this.

References

Batjes, N. H.: ISRIC-WISE global data ser of derived soil properties

on a 0.5 by 0.5 degree grid (version 3.0), ISRIC – World Soil Information, Wageningen, 2005.

Brown, C. D., Dufour‐Tremblay, G., Jameson, R. G., Mamet, S. D., Trant, A. J., Walker, X. J., Boudreau, S., Harper, K. A., Henry, G. H. R., Hermanutz, L., Hofgaard, A., Isaeva, L., Kershaw, G. P., and Johnstone, J. F.: Reproduction as a bottleneck to treeline advance across the circumarctic forest tundra ecotone, Ecography, 42, 137-147, 10.1111/ecog.03733, 2018.

Callaghan, T. V., Jonasson, C., Thierfelder, T., Yang, Z., Hedenas, H., Johansson, M., Molau, U., Van Bogaert, R., Michelsen, A., Olofsson, J., Gwynn-Jones, D., Bokhorst, S., Phoenix, G., Bjerke, J. W., Tommervik, H., Christensen, T. R., Hanna, E., Koller, E. K., and Sloan, V. L.: Ecosystem change and stability over multiple decades in the Swedish subarctic: complex processes and multiple drivers, Philos Trans R Soc Lond B Biol Sci, 368, 20120488, 10.1098/rstb.2012.0488, 2013.

Cleveland, C. C., Townsend, A. R., Schimel, D. S., Fisher, H., Howarth, R. W., Hedin, L. O., Perakis, S. S., Latty, E. F., Von Fischer, J. C., Elseroad, A., and Wasson, M. F.: Global patterns of terrestrial biological nitrogen (N2) fixation in natural ecosystems, Global Biogeochemical Cycles, 13, 623-645, 10.1029/1999gb900014, 1999.

Epstein, H. E., Kaplan, J. O., Lischke, H., and Yu, Q.: Simulating Future Changes in Arctic and Subarctic Vegetation, Computing in Science & Engineering, 9, 12-23, 10.1109/mcse.2007.84, 2007.

Hengl, T., Mendes de Jesus, J., Heuvelink, G. B., Ruiperez Gonzalez, M., Kilibarda, M., Blagotic, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, PLoS One, 12, e0169748, 10.1371/journal.pone.0169748, 2017.

Hofgaard, A., Dalen, L., and Hytteborn, H.: Tree reqruitment above the treeline and potential for climate-driven treeline change, Journal of Vegetation Science, 20, 1133-1144, 2009.

Kollas, C., Vitasse, Y., Randin, C. F., Hoch, G., and Korner, C.: Unrestricted quality of seeds in European broad-leaved tree species growing at the cold boundary of their distribution, Ann Bot, 109, 473-480, 10.1093/aob/mcr299, 2012.

Körner, C., Basler, D., Hoch, G., Kollas, C., Lenz, A., Randin, C. F., Vitasse, Y., and Zimmermann, N. E.: Where, why and how? Explaining the low-temperature range limits of temperate tree species, Journal of Ecology, 104, 1079-1088, 10.1111/1365-2745.12574, 2016.

Rees, W. G., Hofgaard, A., Boudreau, S., Cairns, D. M., Harper, K., Mamet, S., Mathisen, I., Swirad, Z., and Tutubalina, O.: Is subarctic forest advance able to keep pace with climate change?, Glob Chang Biol, 26, 3965-3977, 10.1111/gcb.15113, 2020.

Scherrer, D., Vitasse, Y., Guisan, A., Wohlgemuth, T., Lischke, H., and Gomez Aparicio, L.: Competition and demography rather than dispersal limitation slow down upward shifts of trees’ upper elevation limits in the Alps, Journal of Ecology, 108, 2416-2430, 10.1111/1365-2745.13451, 2020.

Smith, B., Wårlind, D., Arneth, A., Hickler, T., Leadley, P., Siltberg, J., and Zaehle, S.: Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model, Biogeosciences, 11, 2027-2054, 10.5194/bg-11-2027-2014, 2014.

Sturm, M.: Changing snow and shrub conditions affect albedo with global implications, Journal of Geophysical Research, 110, 10.1029/2005jg000013, 2005.

Sturm, M., Holmgren, J., McFadden, J. P., Liston, G. E., Chapin, F. S., and Racine, C. H.: Snow–Shrub Interactions in Arctic Tundra: A Hypothesis with Climatic Implications, Journal of Climate, 14, 336-344, 10.1175/1520-0442(2001)014<0336:Ssiiat>2.0.Co;2, 2001.

Sundqvist, M. K., Björk, R. G., and Molau, U.: Establishment of boreal forest species in alpine dwarf-shrub heath in subarctic Sweden, Plant Ecology & Diversity, 1, 67-75, 10.1080/17550870802273395, 2008.

Truong, C., Palme, A. E., and Felber, F.: Recent invasion of the mountain birch Betula pubescens ssp. tortuosa above the treeline due to climate change: genetic and ecological study in northern Sweden, J Evol Biol, 20, 369-380, 10.1111/j.1420-9101.2006.01190.x, 2007.

Yang, Z., Hanna, E., and Callaghan, T. V.: Modelling surface‐air‐temperature variation over complex terrain around abisko, swedish lapland: uncertainties of measurements and models at different scales, Geografiska Annaler: Series A, Physical Geography, 93, 89-112, 10.1111/j.1468-0459.2011.00005.x, 2011.

Yang, Z., Hanna, E., Callaghan, T. V., and Jonasson, C.: How can meteorological observations and microclimate simulations improve understanding of 1913-2010 climate change around Abisko, Swedish Lapland?, Meteorological Applications, 19, 454-463, 10.1002/met.276, 2012.

Gustafson et al. aim at modelling future treeline position in northern Sweden and the causes that control potential shifts. I am familiar with the associated theory and the region. Using environmental data and biological response functions, a digital vegetation model is applied (with the treeline forming species of this region, Betula pubescence ssp.).

We would like to thank Professor Körner for his comments on our manuscript and for giving us the opportunity to clarify a few aspects of our narrative.

Such a model has a predefined response hierarchy, that is, assumptions on both, the relative importance of drivers and the direction of causality. These assumptions, though absolutely central, are not mentioned upfront, but they become obvious as one reads the text.

All models are simplifications of reality and have assumptions underpinning the processes and parameters included in the model. One of the motivations of this study and the use of such a high-resolution grid was to test some of the assumptions built into the model. The scale used in the study enables comparison to experimental and observational data, which is more difficult to do when, as is typical, the model output represents the average over a large spatial scale. The Abisko region has been thoroughly studied with diverse observational datasets, providing an ideal laboratory for benchmarking a model. A modelling study complements empirical approaches in interpreting the causalities in observed dynamics of treeline advance and vegetation shifts. Lastly, the study aims to simulate future vegetation shifts in the area, something only a model can do.

One of the key assumptions is that these trees are C limited and that photosynthesis, A f (T, PFD, CO2) drives growth. Starting with such assumptions, the inevitable outcome is that CO₂ matters for growth, although it may not matter for treeline position, depending in other assumptions. Yet, in my view this is dressing the horse from the tail. It became obvious in recent decades that outside horticulture and agro-conditions, growth controls A via phloem downloading on demand for C, and this demand is set by meristem activity and other sinks. Not surprisingly manipulating C supply in the field neither rose growth or productivity in alpine vegetation, nor in treeline trees (there were transitory effects on young,  isolated Larix individuals in exceptionally warm summers, that did not affect final biomass data, pine was never affected). None of these works are cited (the only reference to CO₂ experiments is the differential response of two upper montane understorey shrubs).

I do respect the skills of the authors to parameterize and handle such a complex set of algorithms, but the underlaying rationale reflects our understanding of causalities in the 1980s.  I am quite aware that starting with modelling growth rises other issues, but several teams have no engaged (Simone Fatichi, Andrew Friend, several papers) and even the dendro-community is now moving forward in that direction (read e.g. Jan Tumajer et al. Frontiers in Plant Sciences 28 Jan 2021). They are still unable to handle resource supply as modulating factors.

LPJ-GUESS includes representations of several processes that can interact to drive or constrain changes in vegetation composition and ecosystem functioning, and moreover these vary over time in response to the evolving ecosystem state. An example is the treeline and productivity changes investigated in our study. Specifically, the following are tested:

• Climatic – We simulate the influence of historic and future climate (temperature, precipitation and solar radiation) change through a range of scenarios extracted from the CMIP5 project (Taylor et al., 2012).
• Nutrient limitation – This is tested in different nitrogen deposition scenarios (see section 2.4.3 in the manuscript) where we vary the nitrogen load in our projection simulations.
• Productivity – We mainly test the influence of productivity on vegetation shifts through our CO₂-fertilisation experiment (see section 2.4.2 in the manuscript).
• Ecological – A prerequisite for treeline advance is the establishment and growth of trees above the current treeline. Tree seedlings will have to compete with existing vegetation during their establishment (Grau et al., 2012; Lett and Dorrepaal, 2018). Our model includes a representation of interspecific competition for light, soil water and nutrients.

Productivity does not drive treeline advance – we agree with this, in fact it is an outcome of our study, demonstrating that emergent dynamics can not necessarily be trivially predicted from the inputs and modelled processes. Specifically, in our CO₂-only experiment, the treeline did not advance despite increased GPP (Photosynthesis) (See Sec. 3.3.2 & Fig. 5b). If anything, this indicates that the lack of a correlation between tree productivity and treeline advance in the real world could have more than one mechanistic interpretation. There is no scientific consensus yet that meristem activity universally explains treeline position. Indeed it would be strange if evolution were not ‘smart’ enough to make maximal use of available resources to drive fitness.

Our simulation results are not inconsistent with the observation of ample carbon storage in trees close to the treeline (Hoch and Körner, 2012), which indicates that trees close to the treelines do not suffer from carbon shortage. To further highlight this, we have added results for NPP to the results (L266-269; L301-306) and discussion (L404-407) sections of the manuscript. We do acknowledge (Table 1) that the simulated over-estimation of biomass is a limitation of the model and an area where potential improvements can be made. While a few factors (e.g., herbivory from both mammals and insects, mainly Epirrita autumnata) may reconcile the model results with the observations, they cannot fully explain our over-estimation of biomass. Recognising this, and in line with the comments by the reviewer, we put forward temperature limitations on xylogenesis (wood formation) as a potential area for model improvements in the future, as has also been done by others (e.g., Friend et al., 2019; Leuzinger et al., 2013; Pugh et al., 2016) - see our Discussion (see L425-438). This will be of importance for future projections of boreal and Arctic carbon budget estimations and could potentially alter the simulated treeline advance in our Abisko domain. Lastly, our modelled treeline advance is not only constrained by tree physiological factors (e.g., xylogenesis or photosynthesis), but also by ecological factors (e.g., interspecific plant competition and soil nutrient cycling). The importance of including ecological processes in model studies of treeline rather than solely considering bioclimatic limits to treeline advance has also been emphasised by others (e.g., Scherrer et al., 2020). The model does however include hard limits to vegetation distribution through the bioclimatic envelope of each PFT. We have added a few sentences in the methods to clarify this (see L153-156).

We note that more mechanistic treatments of migration and the spread of seeds could also play a role, as could altered disturbance patterns, type and intensity, as was also pointed out by reviewer #1.

Added texts in the manuscript:

Line 153-156:

The bioclimatic envelope is a hard limit to vegetation distribution, intended to represent the physiological niche of a PFT. Furthermore, the climate envelope is a proxy for physiological processes such as meristem activity that may set species ranges, but also for climatic stressors such as tissue freezing.

Line 266-269:

NPP for IBS in the forest increased from 96 gC m^-2 yr^-1 to 180 gC m^-2 yr^-1 over our historic period (1913-2000). Corresponding values at the treeline did not increase but were stagnant at around 60 gC m^-2 yr^-1. Above the treeline, IBS showed very low NPP values (<15 gC m^-2 yr^-1), while NPP for the dominant shrub (LSE) doubled from 20 gC m^-2 yr^-1 to 40 gC m^-2 yr^-1.

Line 301-306:

Forest NPP, wherein IBS was always dominant, increased from 200 gC m^-2 y^-1 in year 2000 to 300 (220 – 375) gC m^-2 yr^-1 and 490 (380 - 610) gC m^-2 yr^-1 for RCP 2.6 and RCP 8.5, respectively, over the projection period. NPP for the same period for IBS at the treeline increased slightly from 60 gC m^-2 yr^-1 to 80 (74 - 90) gC m^-2 yr^-1 and 104 (80 - 116) gC m^-2 yr^-1 for RCP2.6 and RCP8.5. Above the treeline NPP remained low (<25 gC m^-2 yr^-1) for IBS in all scenarios and always had a lower NPP than the most productive shrub PFT (LSE). NPP for this shrub was 49 (24 - 64) gC m^-2 yr^-1 and 130 (81 – 180) gC m^-2 yr^-1.

Line 404-407:

Furthermore, although NPP for IBS was lower at the treeline than in the forest, it was never close to zero. Such a pattern, which was seen above the treeline, would have indicated a stagnant growth, and that the productivity and carbon costs of maintaining a larger biomass would have cancelled each other out.

Line 441 – 454:

IBS at the treeline had a positive carbon balance (NPP) and was thus not directly limited by its productivity in our simulations. This is consistent with observations of ample carbon storage in treeline trees globally (Hoch and Körner, 2012). The modelled treeline is thus not set by productivity directly but rather by competition, as other PFTs become more productive above the treeline. Whether the treeline is set by productivity constraints or by cold temperature limits on wood formation and meristematic activity has been a subject of some discussion (Körner, 2015, 2003; Körner et al., 2016; Fatichi et al., 2019; Pugh et al., 2016). DVMs have traditionally assumed photosynthesis to be constraining for growth, and thus species distributions. On the other hand, trees close to the treeline have not shown any shortage of carbon for growth (Hoch and Körner, 2012). Furthermore, enhancement of photosynthesis through added CO₂ has also not always resulted in increased tree growth close to the treelines (Dawes et al., 2013), and wood formation is slow around 5°C, leading to a hypothesis of reversed control of plant productivity and range distributions of trees (Körner, 2015). Lately, ecological interactions as a component in the control of treeline position, rather than just considering hard limits to species distributions, has been a subject of more attention in modelling studies (See for ex., Scherrer et al., 2020). These studies add an extra dimension to the discussion as they do not only consider plant physiology but also broadly accepted ecological concepts such as realised versus fundamental niches.

Line 457-465:

… elevation than observed by Hedenås et al. (2011). The mean annual biomass carbon increase in the same dataset is, although highly variable, on average 2.5 gC m^-2 yr^-1 between 1997 and 2010. As simulated GPP and LAI were within the range of observations in the area (Rundqvist et al., 2011; Ovhed and Holmgren, 1996; Olsson et al., 2017), this indicates a coupling between photosynthesis and growth in the model that is stronger than observed. Terrestrial biosphere models often overestimate biomass in high latitudes (Pugh et al., 2016; Leuzinger et al., 2013) and potentially lack processes that likely limit growth close to low temperature boundaries. Examples of such processes are carbon costs of nitrogen acquisition (Shi et al., 2016), including costs for mycorrhizal interactions (Vowles et al., 2018), and temperature limits on growth increment (Friend et al., 2019), i.e., decoupling of growth and photosynthesis. However, data on carbon allocation and its temperature dependence is scarce (Fatichi et al., 2019).

With these concerns, the results of the modelling reflects the assumptions. If one assumes soil fertility matters for treeline trees and selects N to represent these nutrients, the outcome is that N matters. 

Dynamic vegetation models (DVMs) are designed to test multiple and interactive driving factors and processes under different environmental conditions. LPJ-GUESS (Smith et al., 2001; Smith et al., 2014), does not build on the assumption of any single driving or modulating factor. The importance of each driver is an emergent outcome of the simulated dynamics in response to variation in the drivers, and the evolving state of the system (e.g. soil N availability, plant community structure) which also modify the processes. Thus we respectfully disagree with the statement “the results of the modelling reflects the assumptions“. While this is true by definition for any model, the point the reviewer is making, that the simulated treeline dynamics could be trivially pre-empted, knowing the input to the model and the process formulations used, is not correct.

If soil fertility were controlling treeline position there should not be a global treeline isotherm and treeline should be at higher elevation on good soils and at lower elevations on poor soils, not what we see in the field (e.g. soils developing on young glacial deposits versus treelines on geologically old, weathered, low latitude mountains.

The reviewer’s comment states that, given our conclusions, local variations in soil fertility would give rise to differences in treeline elevation and that no such pattern is seen in global treeline records. This would be true if no recognition of temperature as a controlling factor were made. We do not dispute the importance of temperature as an overall control on treeline position, but stand by our conclusion, consistent with physiological, demographic, and ecological principles and assumptions around which there is broad scientific consensus, that nitrogen limitation will constrain rates of future treeline advance in this area.

The global correlation between treelines and the 6.7°C isotherm (6-7°C in arctic/boreal regions; Körner and Paulsen, 2004) is not disputed in this study. In fact, results from our historic simulations corroborate this, at least for this region since the position of the simulated treeline correlates well with an isotherm close to this limit (Fig S1.4; supplementary materials). The lag arises in our future projections of treeline advance when the isotherm displacement is more rapid than the treeline advance in most scenarios (Fig 4). We traced this effect to N-limitation in the soils as we see a faster treeline migration in our climate change only simulations compared to the future projections when N-deposition is decreasing.

Experiments by Hoch (2013) revealed that there is no compensatory effect of nutrient addition to low temperature constraints of growth. There is also no direct link between tree vigor and treeline position. Trees at treelines in Bolivia and Tibet at close to 5000 m elevation hardly grow (minute tree ring width), because they are clearly moisture and thus, nutrient limited.

The model simulates induced treeline advance when only nitrogen deposition increased (Fig 5c). This was a result of an increased allocation to above-ground biomass with a subsequent advantage in the light competition. The reviewer refers to the study by Hoch (2012) as evidence against such an effect. However, the Hoch study was performed in a fully factorial design with seedlings grown in phytotrons with controlled environment at 6 ^oC or 12 ^oC and half of the seedlings fertilised. The cold temperature seedlings showed no response to added fertiliser, while warm temperature seedlings showed a strong biomass increase with fertilisation. We would argue that such a controlled climate study may not mimic the role of nitrogen effects in a full ecosystem setting.

I hope these comments are useful for revisiting the rationales underpinned in this model. I read the other report in copernicus, thanks for providing it. It seems to address additional facets of treeline formation, but does not touch upon the more fundamental bias regarding the assumptions that drive the model output. 

We would like again to thank Professor Körner for his comments. We would also argue that our study makes a meaningful contribution to the treeline literature, and hope that our findings will stimulate further experimental and modelling studies of treeline advance and associated feedbacks to climate.

References

Dawes, M. A., Hagedorn, F., Handa, I. T., Streit, K., Ekblad, A., Rixen, C., Körner, C., and Hättenschwiler, S.: An alpine treeline in a carbon dioxide-rich world: synthesis of a nine-year free-air carbon dioxide enrichment study, Oecologia, 171, 623-637, 10.1007/s00442-012-2576-5, 2013.

Fatichi, S., Pappas, C., Zscheischler, J., and Leuzinger, S.: Modelling carbon sources and sinks in terrestrial vegetation, New Phytol, 221, 652-668, 10.1111/nph.15451, 2019.

Friend, A. D., Eckes-Shephard, A. H., Fonti, P., Rademacher, T. T., Rathgeber, C. B. K., Richardson, A. D., and Turton, R. H.: On the need to consider wood formation processes in global vegetation models and a suggested approach, Annals of Forest Science, 76, 10.1007/s13595-019-0819-x, 2019.

Grau, O., Ninot, J. M., Blanco-Moreno, J. M., van Logtestijn, R. S. P., Cornelissen, J. H. C., and Callaghan, T. V.: Shrub-tree interactions and environmental changes drive treeline dynamics in the Subarctic, Oikos, 121, 1680-1690, 10.1111/j.1600-0706.2011.20032.x, 2012.

Hedenås, H., Olsson, H., Jonasson, C., Bergstedt, J., Dahlberg, U., and Callaghan, T. V.: Changes in Tree Growth, Biomass and Vegetation Over a 13-Year Period in the Swedish Sub-Arctic, Ambio, 40, 672-682, 10.1007/s13280-011-0173-1, 2011.

Hoch, G.: Reciprocal root-shoot cooling and soil fertilization effects on the seasonal growth of two treeline conifer species, Plant Ecology & Diversity, 6, 21-30, 10.1080/17550874.2011.643324, 2012.

Hoch, G. and Körner, C.: Global patterns of mobile carbon stores in trees at the high-elevation tree line, Global Ecology and Biogeography, 21, 861-871, 10.1111/j.1466-8238.2011.00731.x, 2012.

Körner, C.: Carbon limitation in trees, Journal of Ecology, 91, 4-17, 2003.

Körner, C.: Paradigm shift in plant growth control, Curr Opin Plant Biol, 25, 107-114, 10.1016/j.pbi.2015.05.003, 2015.

Körner, C. and Paulsen, J.: A World-Wide Study of High Altitude Treeline Temperatures, Journal of Biogeography, 31, 713-732, 2004.

Körner, C., Basler, D., Hoch, G., Kollas, C., Lenz, A., Randin, C. F., Vitasse, Y., and Zimmermann, N. E.: Where, why and how? Explaining the low-temperature range limits of temperate tree species, Journal of Ecology, 104, 1079-1088, 10.1111/1365-2745.12574, 2016.

Lett, S. and Dorrepaal, E.: Global drivers of tree seedling establishment at alpine treelines in a changing climate, Functional Ecology, 32, 1666-1680, 10.1111/1365-2435.13137, 2018.

Leuzinger, S., Manusch, C., Bugmann, H., and Wolf, A.: A sink-limited growth model improves biomass estimation along boreal and alpine tree lines, Global Ecology and Biogeography, 22, 924-932, 10.1111/geb.12047, 2013.

Olsson, P.-O., Heliasz, M., Jin, H., and Eklundh, L.: Mapping the reduction in gross primary productivity in subarctic birch forests due to insect outbreaks, Biogeosciences, 14, 1703-1719, 10.5194/bg-14-1703-2017, 2017.

Ovhed, M. and Holmgren, B.: Modelling and measuring evapotranspiration in a mountain birch forest, Ecological Bulletins, 45, 31-44, 1996.

Pugh, T. A. M., Muller, C., Arneth, A., Haverd, V., and Smith, B.: Key knowledge and data gaps in modelling the influence of CO2 concentration on the terrestrial carbon sink, J Plant Physiol, 203, 3-15, 10.1016/j.jplph.2016.05.001, 2016.

Rundqvist, S., Hedenås, H., Sandström, A., Emanuelsson, U., Eriksson, H., Jonasson, C., and Callaghan, T. V.: Tree and Shrub Expansion Over the Past 34 Years at the Tree-Line Near Abisko, Sweden, Ambio, 40, 683-692, 10.1007/s13280-011-0174-0, 2011.

Scherrer, D., Vitasse, Y., Guisan, A., Wohlgemuth, T., Lischke, H., and Gomez Aparicio, L.: Competition and demography rather than dispersal limitation slow down upward shifts of trees’ upper elevation limits in the Alps, Journal of Ecology, 108, 2416-2430, 10.1111/1365-2745.13451, 2020.

Shi, M., Fisher, J. B., Brzostek, E. R., and Phillips, R. P.: Carbon cost of plant nitrogen acquisition: global carbon cycle impact from an improved plant nitrogen cycle in the Community Land Model, Glob Chang Biol, 22, 1299-1314, 10.1111/gcb.13131, 2016.

Smith, B., Prentice, I. C., and Sykes, M. T.: Representation of vegetation dynamics in the modelling of terrestrial ecosystems: comparing two contrasting approaches within European climate space, Global Ecology and Biogeography, 10, 621-637, 10.1046/j.1466-822X.2001.t01-1-00256.x, 2001.

Smith, B., Wårlind, D., Arneth, A., Hickler, T., Leadley, P., Siltberg, J., and Zaehle, S.: Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model, Biogeosciences, 11, 2027-2054, 10.5194/bg-11-2027-2014, 2014.

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, Bulletin of the American Meteorological Society, 93, 485-498, 10.1175/bams-d-11-00094.1, 2012.

Vowles, T., Lindwall, F., Ekblad, A., Bahram, M., Furneaux, B. R., Ryberg, M., and Bjork, R. G.: Complex effects of mammalian grazing on extramatrical mycelial biomass in the Scandes forest-tundra ecotone, Ecol Evol, 8, 1019-1030, 10.1002/ece3.3657, 2018.