Manuscript: gmd-2021-422

Response to Referees

Dear Dr. Sato,

Thank you for giving us the opportunity to address the comments to the manuscript “Tree migration in the dynamic, global vegetation model LPJ-GM 1.0: Efficient uncertainty assessment and improved dispersal kernels”. We would like to thank you and the reviewers for taking the time to provide feedbacks on our manuscript. We have been able to incorporate most of the suggestions made by the reviewers. Please see below the point-by-point response to the reviewers’ comments. All line numbers in the Response refer to the revised manuscript.

Anonymous Referee #1

Major Comment 1: The study appears to focus on migration of European tree species, but this geographical focus is not mentioned very explicitly. E.g., in the review of literature and observational data (2.1), there appears to be a focus on Europe, but this is not presented until the methods section. I would recommend to mention the European focus early on, also in the abstract and probably in the title.

Response: Thank you for your suggestion. We added explicit references to European forests in the title (“Tree migration in the dynamic, global vegetation model LPJ-GM 1.1: Efficient uncertainty assessment and improved dispersal kernels of European trees”), the Abstract (L20): “This approach was tested on the newly-implemented migration module of the state-of-the-art DGVM LPJ-GM, focusing on European forests.”, and in the Introduction (L120): “In this study, we present an efficient uncertainty assessment of model selection and parameter estimates for the newly-implemented migration module of LPJ-GM. We focus our assessment on European forests as estimates of migration rates from paleo-records are available for most major European tree taxa (see Sect. 2.1).”

Major Comment 2: The choice of the simulation setup is not described in great detail, and is not justified in the text. Does the resolution in the simulation (pixels of 1x1 km) fit with the dispersal distances? Intuitively, my feeling would be that this would require a higher resolution to be resolved, in particular regarding the substantial differences in the tested kernels in the range 0-400 m (Fig. 3), and also in the light of the expected migration rates (Table 1) and the resulting dispersal distances (Fig. 6). I welcome the authors to comment on that and justify their choice of the resolution in the methods section. Apart from that, the authors should explain how the “static suitable climate” (l. 254) was determined. Is this the same climate for all species?

Response: Thank you for highlighting these points for further clarification. We added more details in the section where the simulation setup is described (Sect. 2.3 Simulation protocol).

Concerning the explanation for the climatic input, we added the following paragraph (L301-304): “We selected a spatial and temporal subset from the climatic dataset provided by Armstrong et al. (2019) so that temperature was not limiting the survival and establishment for all species according to the species-specific bioclimatic limits defined by the LPJ instruction file (i.e. the minimum month mean temperature for survival, tcmin_surv; see Code and data availability and Supplement 5). We applied this suitable climate to all species and an entirely permeable terrain to all grid cells and across all simulation years in order to reproduce optimal environmental conditions for tree migration.”

Concerning the resolution of the simulation grid cell and the dispersal patch size, we added details on the spatial simulation domain, incorporating suggestions from Comment 9 by the Friendly Reviewer listed below (L286-297): “The climate grid of 0.5 degrees resolution was subdivided into smaller cells of 1 km² area, where vegetation dynamics are simulated at a patch level of 1000 m², as is usually done for LPJ-GUESS simulations. Simulations were performed for a total of 500 years, covering an area of 201 x 201 cells with corridors located on the perimeter and the two major diagonals of the domain for a total of 1,197 simulated cells (see Fig. S1). After the spin-up phase, migrating tree species were allowed to establish freely in the upper-left corner of the simulated landscape (the starting point of migration). In the dispersal routine, seeds were dispersed at the beginning of each year according to the kernel function (Eq. (A3)), which was applied from the centres of each 1 km² cell (thus, the minimum dispersal distance  to reach a different sink cell is 500 m starting from the middle of the source cell). The distance between neighbouring cells is surpassed by the long distance seed dispersal and, though the average distance of local seed dispersal (SDD) is below this distance, the overall dispersal kernel still stretches outside of the source cell with a SDD of 25–100 m (see Table B3), thus allowing migration to neighbouring cells. To additionally reduce the error introduced by discretizing the seed dispersal kernel to 1 km², the kernel was first computed on the finer resolution of 100 x 100 m cells and subsequently summed up over the 1 km² cells.”. Additionally, we added a sentence in the Result section with reference to the above paragraph (L551-552): “Regarding the influence of the dispersal patch size, our parametrization holds for the spatial resolution of 1 km (with a finer resolution of 100 m for the simulation of the dispersal kernel) (see Sect. 2.3).”. Thus, the actual resolution of the simulation of seed dispersal is 100 x 100 m among the finer cells, which we consider sufficient to compute the seed dispersal.

Major Comment 3: And finally a question about the numbering of model versions: In my understanding, the original publication (Lehsten et al. 2019) presents version 1.0, and the current manuscript describes version 1.1 (l. 129). For clarity, I would recommend to use the new version number in the title, rather than referring to it as an update of version 1.0.

Response: Thank you for mentioning this issue. Indeed, nomenclature is very important for the clusters of LPJ-models so that it will be advisable to distinguish our current model from the previously published one. We have now modified the title following your suggestion.

Minor Comment 1: l. 19 remove comma after “DGVM”.

Response: The comma has now been removed (L20).

Minor Comment 2: l. 80 “set” is singular, I recommend to alter “minimize” to “minimizes”.

Response: We modified “minimize” to “minimizes” (L89).

Minor Comment 3: l. 97 replace “little” by “small”.

Response: We replaced “little” by “small” (L106).

Minor Comment 4: l. 98 explain what “seed limitation” means in a dispersal context.

Response: We added the following explanation for seed limitation (L108-109): “…with and without dispersal limitation, i.e. seed distribution was simulated by dispersal from mature trees, or seeds were distributed freely across the whole simulation domain, respectively”.

Minor Comment 5: l. 163-168 Are examples 1-4 all originating from MacDonald (1993)? This is not entirely clear.

Response: We updated this section with references and following the suggestion from Comment 3 by the Friendly Reviewer listed below (L172-180).

Minor Comment 6: l .196 replace “above-mentioned” by “mentioned above”.

Response: The sentence was replaced following the suggestion from Comment 7 by the Friendly Reviewer listed below (L172-180).

Minor Comment 7: explain the computation of SI – Is delta y computed as difference of all combinations of 2 points (out of the 5 points sampled)? And how does this match with the 25% value change of delta x? Or is it only the difference between subsequent points? In case of the latter, it would be enough to compute the difference between the outcomes of minimum and maximum of the distribution, divided by the full range of the distribution.

Response: Thank you for pointing this out. The deltas (x and y) of SI are the differences between subsequent points. Following your suggestion, the formulation of SI was modified as (y_max – y_min)/(x_imax – x_imin) (L346-351). The re-formulation does not affect the values of SI.

Minor Comment 8: Eq. 5 Equation appears to lack a summation in order to compute residuals. Also, it is unclear why only integer values (l. 329) are used here – migration rate is a continuous variable, right? In more general terms, the use of all integer values between lower and upper boundaries appears to assume a uniform distribution between these extremes. I would suggest that this assumption is mentioned explicitly.

Response: Following suggestions from Comments 7, 12, 14 and 20 by the Friendly Reviewer listed below, we changed the formula of species-specific residuals by using only maxima of migration rate estimates as observed values (L380-383). With this approach, there is no use of integers (i.e. values between minimum and maximum estimates of migration rates). There is no summation in this formula (if compared to Eq. (7)) as residuals are species-specific and only one value of maximum migration rate is available per species.

Minor Comment 9: l. 345 are the fat-tailed kernels implemented to represent the long-distance dispersal only (so only the second pdf in Eq. A3), or do you use them to represent both? And are the results displayed in Fig. 3 and Fig. 6 representing the LDD part only, or the combination of both?

Response: Fat-tailed kernels are used for both short-distance and long-distance dispersal. We added this clarification in the relevant parts of the text as you suggested: 1) in Sect 2.6.1 (L407-410): “we… ran simulations with default parameter values for the two linearly combined pdfs as calculated from the average dispersal distances at the local scale (SDD_d) and for long-distance dispersal events (LDD_d) (i.e. fat-tailed functions are implemented to represent the pdfs of both SDD and LDD components in Eq. (A3);…)”; 2) in the explanation of Eq. (A3) (L789-790): “  and  are the scale parameters calculated from  and , respectively (the average distances for SDD and LDD; see Table 3)”; in the legend of Fig. 3 (L559): “Kernel formulae are implemented as pdfs in Eq. (A3) for both the SDD and LDD component.” and Fig. 6 (L574): “Each species-specific kernel function is implemented as pdf in Eq. (A3) for both the SDD and LDD component.”.

Minor Comment 10: l. 353 ensure that you mention earlier that you have selected five kernel functions – it is only clear from the table which five functions are meant here.

Response: We added a reference to the five functions before the table, when fat-tailed kernels are mentioned for the first time (L399-401): “we decided to implement five additional dispersal kernels to better represent long-distance dispersal (LDD) events in the migration model: exponential power (ExpPow), Weibull, bivariate Student’s t (twoDt), logistic, log-hyperbolic secant (LogSec) (Table 3).”.

Minor Comment 11: Eq. 6 Add a reference for the equation, or explain the equation in more detail (e.g., where the square root of 2 comes from).

Response: As suggested, we added the following sentence to clarify Eq. (6) (L436-438): “Eq. (6) was derived from Eq. (3) (Snowling and Kramer, 2001) by setting the weighting constants for sensitivity and error to 1 (i.e. error and sensitivity were valued equally relative to each other in the calculation of the model utility).”.

Minor Comment 12: l. 391 replace “variable” with “variables”.

Response: We replaced “variable” with “variables” (L449).

Minor Comment 13: Eq. 8 and its explanation (l. 391) ensure that subscripts are corrected – p and r appear to be functions of both i and l.

Response: Thank you for spotting the error. Eq. (8) and its explanation (L448-450) are now modified following your suggestion and the original formulation in Snowling and Kramer (2001), Eq. (1).

Minor Comment 14: l. 411 replace “lead” with “leads”.

Response: We replaced “lead” with “leads” (L469).

Minor Comment 15: l. 574 replace “high leptokurtosis” with “high kurtosis” or “leptokurtic”.

Response: We replaced “high leptokurtosis” with “high kurtosis” (L643).

Minor Comment 16: l. 590 replace “fat-tail” with “fat-tailed” for consistency.

Response: We replaced “fat-tail” with “fat-tailed” (L659).

Anonymous Referee #2

Comment 1: My only question for the authors to consider has to do with how they address the role of anthropogenic landcover fragmentation in vegetation range shifts. The authors go into substantive detail about reproducing migration rates since the last glacial maximum, and explore the essential ingredients of their model necessary to reproduce post-glacial migration rates. But – and the authors acknowledge this – post-glacial migration rates manifested in a landscape not yet fragmented by human land-use and built environments. So while the potential effects of habitat fragmentation (~L145) and the potential limitations of this model for predicting future range shifts (~L665) get some mention, anthropogenic forcing struck me as under-emphasised, overall. In terms of predicting future range shifts at continental/regional scales, fragmentation and built environments function like natural physical barriers (e.g., mountain ranges), and looking forward, might reasonably be an external forcing as strong as climate change. Do post-glacial vegetation migration models work in the latest Holocene? Will understanding range shifts through aggressively fragmented landscapes require a different kind of model?

Response: Thank you for mentioning this important aspect of range shift dynamics. Indeed, as you mentioned, it is predicted (and it has been observed since decades) that anthropogenic effects (either as landscape alteration or assisted migration) will have a significant impact on species’ range shift. As anthropogenic effects are usually challenging to disentangle from other drivers, we first decided to explore the more natural (i.e. ecological) processes underlying plant migration by simulating vegetation expansion under optimal conditions (with optimal climate and no barriers, including anthropogenic disturbances). We are currently simulating the forest expansion of Europe throughout the whole Holocene (starting from 20,000 ya) by including migration with the optimal parameters we obtained from the current manuscript. Simulations would then be validated by pollen-data and one of our aims is to understand whether the performance of LPJ-GM 1.1 is higher in prehistoric times with respect to Late Holocene when human disturbance starts to increase. This might give an insight on the anthropogenic effect on vegetation expansions as in the current simulations we include no anthropogenic effect. However, the LPJ-model family has an option to include anthropogenic effect by using different land use covers (and thus fragmentation and habitat loss by urban and agricultural areas) and terrain permeability in the simulation of migration. We are also planning to run future simulations of vegetation range shifts in Europe by including land use change to test for the anthropogenic effect on plant migration by the end of the 21^st century (this is now briefly mentioned in the manuscript as a future avenue; see below).

Comment 2: I don't suggest the authors change anything about the quantitative exercise they have undertaken. My only query is whether they might amend a portion of their introduction to push the problem of anthropogenic fragmentation, and the open questions that raises for global vegetation migration (real and modelled), more into the foreground. Doing so would shift the framing of anthropogenic disturbance that comes at the end of the Discussion away from being a "limitation" and toward an avenue of future research.

Response: Following your suggestions, we added a paragraph early in the Introduction (L36-47) concerning the effect of anthropogenic disturbances on plant migration: “ Nevertheless, contemporary range shifts are expected to differ from past dynamics, as species are submitted to different conditions, such as the higher velocity of current global warming relative to past (post-glacial) climate changes and the limited anthropogenic influence on the landscape in prehistoric times (Nogués-Bravo et al., 2018). Concerning the human impact on landscape reconfiguration, habitat loss has been shown to have a significant and negative impact on vegetation expansion both in models (Collingham and Huntley, 2000; Dullinger et al., 2015; Saltré et al., 2015) and real case studies in the current century (Guo et al., 2018). […]. Additionally, analyses of past vegetation dynamics should aim to identify relevant processes underlying plant range dynamics that can be still valid for future scenarios with increasing anthropogenic alterations of the landscape (Nogués-Bravo et al., 2018).”.

We also mentioned a future avenue for research that includes the anthropogenic effect on plant range shifts in the section of Limitation and Future Challenges (L712-723): “As an avenue of future research, settings with heterogeneous terrains and climate can be applied to simulations to assess the relative contribution of climate, landscape heterogeneity, habitat loss and dispersal ability in the context of real world scenarios. […]. For a more explicit assessment of landscape configuration on species’ range shifts, future simulations can be performed by driving LPJ-GM with projections of climate and land use change throughout the whole 21st century as human-driven habitat loss has already been reported to have significantly affected plant migration at the start of the current century (Guo et al., 2018). In this regard, LPJ-GM can optionally simulate heterogeneous landscapes by using a spatially explicit seed dispersal permeability value (Lehsten et al., 2019), where landscape permeability can be informed by future projections of land use change.”.

Friendly Reviewer #1

Comment 1: L12. It would be preferable, and more accurate, to use 'propagule' rather than seed, as many of the dispersed entities are fruits (e.g. in the case of Betula).

Response: Thank you for pointing out this inaccuracy. We added ‘propagule’ in parenthesis after seed in the Abstract (L14). We decided to keep ‘seed’ as we use this term in the explanation of the model (Method; Appendix).

Comment 2: L147-160: “the role of climate in controlling the rates of tree expansion after the LGM is generally assumed to be less influential than intrinsic dispersal ability as indicated by the wide-spread phenomenon of migration lag (i.e. the delayed arrival of a species into a newly-suitable habitat; Feurdean et al., 2013; Huntley et al. 2013; Giesecke and Brewer, 2018), thus allowing to reduce the number of influential migration drivers during simulation.”. Migration lagging behind climate is only an issue when climatic changes are extremely rapid, and especially when the direction of change is soon afterwards reversed. This was the case during the last glacial stage when Heinrich Events, especially, but also Dansgaard-Oeschger fluctuations, resulted in very rapid climatic changes, but to new conditions that were then short-lived. Once the last termination was passed, and the Holocene was established, climatic changes were orders of magnitude less rapid, and the 'new' conditions were persistent.  As a result, whilst some migration lag may have occurred in the first few millennia of the Holocene, and even that is debatable, thereafter the evidence is pretty clear that climate was the major driver of the continuing migrations and range shifts of tree taxa in Europe, in eastern North America, etc. (see e.g. Prentice et al., 1991).

Prentice, I.C., Bartlein, P.J. & Webb, T., III (1991) Vegetation and climate change in Eastern North America since the last glacial maximum. Ecology, 72, 2038-56.

Response: For the analyses of this paper, we do not assume that the relative influence of dispersal limitation was greater than the climate in determining the range shift of tree species after the LGM. We decided to simulate tree migration in a non-limiting climate so to reduce as much as possible the influence of the environment on the simulated migration rates. Thus, we could assess the impact of the migration parameters and the model structure on the simulated migration speed and possibly optimize the parameters to match the maximum values from the literature (see also response to Comment 22). Following your suggestion, we eliminated the sentence above.

Comment 3: L163-164: “1) Radiocarbon dating of fossil pollen generally provides a coarse temporal resolution, ranging from several decades to more than a century.”. This is essentially irrelevant to the migration rate estimates made from palynological data, because typically the rates are computed by examining isopoll displacements over intervals of a millennium, or more, or using isochrone displacements, where the isochrones are drawn typically at millennial intervals.  The uncertainty in radiocarbon dating for the Holocene, which is the interval for which the migration rate estimates are made, is rarely more then a century or two, and certainly makes a much smaller contribution to the uncertainty in location of isopolls or isochrones than does the relatively sparse spatial distribution of sites with pollen diagrams, or indeed the relatively temporally sparse sampling of many sites with pollen diagrams.

Response: Thank you for making this clarification. We modified the paragraph accordingly (L174-180): “Similarly, radiocarbon dating of fossil pollen generally provides a coarse temporal resolution, ranging from several decades up to more than a century in the post-LGM period (MacDonald, 1993). However, migration rates after the LGM are usually estimated over intervals of a millennium (e.g. Huntley and Birks, 1983; Giesecke et al., 2017), thus the relatively coarse temporal resolution in the radiocarbon dating may not contribute much to the uncertainty in migration estimates. More relevant sources of uncertainty are the uneven spatial distribution of sites with available fossil pollen and the temporally limited sampling of some of these sites (Huntley and Birks, 1983; B. Huntley, 2022, pers. comm., January).”.

Comment 4: L174-176: “However, previous studies based on plant macrofossils and potential glacial tree distribution (e.g. Stewart and Lister, 2001) have also hypothesized the presence of northern refugia during the LGM (above 45º N), which would yield lower rates of northward tree migration (Feurdean et al., 2013)”. These so-called 'northern' or 'cryptic' refugia have been hypothesised, but to my knowledge no convincing evidence has been presented for the distribution of temperate tree species north of 45N in western Europe during the glacial maximum, and even in central and eastern Europe the evidence is far from solid. One serious problem with much of the evidence is with the dating of material, and also a strong tendency for authors to assume presence in areas where a tree was present many millennia before the LGM and then re-appeared in the early Holocene (see Huntley, 2014).  A second problem is the failure of many authors to distinguish between boreal and temperate trees, just lumping all trees together when discussing putative northern refugia.

Huntley, B. (2014) Extreme temporal interpolation of sparse data is not a sufficient basis to substantiate a claim to have uncovered Pleistocene forest microrefugia. The New phytologist, 204, 447-9.

Response: Thank you for proving information and references to the question of the use of cryptic refugia in the estimation of migration rates. We added this information at the end of the section (L222-226): “Finally, we decided to exclude estimates calculated assuming the presence of northern glacial refugia (Feurdean et al., 2013) based on the unlikely survival of temperate tree taxa north of 45º N during the LGM (Tzedakis et al. 2013) and on the lack of strong evidence for the presence of northern refugia (e.g. detection based on extreme temporal interpolation of sparse data or due to the overestimation of occurrences with presence-only spatial distribution models; Huntley, 2014).”

Comment 5: L183: “…established individual, which can be identified at the species level (Binney et al., 2009).”. Although some tree macrofossils can be identified to species level, most are limited to species groups, sub-genera or even entire genera.  Notably, I do not recollect ever having seen Quercus spp. identified from wood, leaves or fruits, except in the British Isles where it is assumed that the only species to be considered are Q. robur and Q. petraea.  Similarly, the wood of the various species of Ulmus cannot be distinguished, and the fruits have rarely if ever been found as macrofossils.  Wood of Salix and Populus is generally considered indistinguishable, and the species within these genera could only potentially be distinguished if very well preserved leaf macrofossils were found, and even then I would consider any identification to species level on that basis as extremely uncertain.  The fruits of the two genera rarely if ever are preserved.  Bud scales just might be distinguishable, but are rarely found or identified to any fine taxonomic level. Even the charcoal upon which Binney places so much evidence has serious limitations in terms of species-level identification.  Pinus charcoal cannot be identified to any finer taxonomic level than can Pinus pollen, i.e. haploxylon vs. diploxylon.

Response: Thank you again for this in-depth information. We modified the text accordingly by deleting “which can be identified at the species level”.

Comment 6: L187-189: “Upper and lower boundaries for the value ranges of migration rates were derived from different empirical studies based on the method employed for their estimation. Pollen-based estimates of maximum rates of spread for common European tree taxa were first summarized by Huntley and Birks (1983). Giesecke et al. (2017) revisited these estimates by correcting for the uneven distributions of pollen diagrams using interpolated maps of pollen percentages and threshold values to reduce the risk of false presence.”. As Huntley & Birks (1983) made clear in their text, they applied threshold pollen values that were assessed as likely to indicate local presence of the taxon, as opposed to long-distance dispersal of pollen.  These threshold values were determined on the basis of mapping of pollen values in surface samples (maps are presented in their book) and a knowledge of the actual distribution of the taxon.  As written, your text implies that Huntley & Birks (1983) did not use threshold values, and that threshold values were first used in the work of Giesecke et al. (2017).  This is clearly not the case.

Response: Thank you for highlighting this error. The text was modified as follows (L200-205): “Pollen-based estimates of maximum rates of spread for common European tree taxa were first summarized by Huntley and Birks (1983) using high abundance thresholds to distinguish the local presence of a taxon from pollen depositions due to rare long-distance dispersal (LDD) events, thus reducing the risk of false presence. Giesecke et al. (2017) revisited these estimates by using over 780 pollen diagrams stored in the European Pollen Database (EPD) along with interpolated maps of pollen percentages and threshold values to infer the distribution and abundances of major European tree taxa in the last 15,000 years.”.

Comment 7: L195-201: “This criterion was based on the assumption that pollen analysis, as above-mentioned, may overestimate the true rates of migration (i.e. determined by first arrival followed by population expansion), especially in the case of tree taxa with long-distance pollen dispersal (LDD). That is, the detection of pollen on the spreading front might indicate an event of first arrival with tree establishment (from which a true rate of migration is calculated) as well as an LDD event with no further establishment and/or population expansion (MacDonald, 1993; Giesecke and Brewer, 2018). In the latter scenario, pollen analysis will 200 result in an apparent rate of migration, which will be an overestimation of the true rate.”. I would argue that it is equally, if not more, likely that (conventional) pollen analysis will not detect the first arrival of a taxon in a region, but rather that pollen abundance values will respond by increasing as the population of the tree expands in the region.  As a result, even the most rapid rates of migration assessed from pollen analytical data are at least as likely to be UNDER-estimates of the true migration rate, if that 'true' rate is based upon the first arrival of the taxon in a region. On this basis, I would argue that any model of dispersal/migration must simulate a rate at least comparable to the maximum rate inferred from pollen analytical data.  Slower rates observed in the palaeoecological record will indicate some 'external' limitation of this rate, most usually by a relatively slower rate of change of climate to the conditions favourable for the species, whether directly (i.e. related to the taxon's physiological requirements/tolerances) or indirectly (by modulating the outcome of some biotic interaction, be it with a competitor, a pathogen, or whatever). It is important to realise, when looking at rates such as those assessed by Huntley & Birks (1983), that threshold pollen abundance values have been used that are selected so as to minimise the likelihood of inferring presence when in reality the pollen is derived from long-distance dispersal, and that the pollen data used are derived from sites where the pollen catchment is at least dominantly extra-local, and is usually regional.  In order to detect first presence of a taxon in an area one would need to examine (a very large number of) small basins where the pollen is derived mostly from the immediate surroundings of the site.  The influence that this can have upon the assessment of a date of 'arrival' is well illustrated by our recent study of Pinus sylvestris forest establishment in Scotland during the Holocene (Huntley & Allen, 2021), in which the date of forest establishment (i.e. a rapid increase in the local population, as reflected by a rapid increase in pollen abundance) assessed from individual small basins can differ from that inferred from the record in a larger basin in the same landscape by as much as half a millennium or more.

Huntley, B. & Allen, J.R.M. (2021) Holocene expansion of the Caledonian pinewoods: spatial and temporal patterns at regional and landscape scales. Plant Ecology & Diversity, 1-24.

Response: Thank you again for clarifying this very important point. Following your comment, we decided to run further analyses (Evaluation of dispersal kernel and Error and Uncertainty Analysis; Sect. 2.6) by considering only the maxima of migration rate estimates as observations. As these new analyses had overall similar results (though not always) to our previous one (which considered either the full range of estimates or the 75^th percentile), we decided to consider only results generated with maxima for the model evaluation (and optimization) of LPJ-GM.

Following the suggestions of this comment, we also modified the text as follows (L214-222): “We then assumed that these ranges might reflect estimates made at different points in space and time of a species’ range during expansion, rather than the uncertainty linked to migration rate estimates (B. Huntley, 2022, pers. comm., January). Thus, the upper boundaries based on pollen analysis should represent the maximum rates achieved by a species during its post-glacial expansion, likely under non-limiting environmental conditions. Additionally, we consider that pollen analyses using abundance thresholds (e.g. Huntley and Birks, 1983) may not detect the first arrival of a taxon at a site and its subsequent establishment as long as it stays regionally rare, but rather the population expansion at a site whether the taxon is previously present in the region or not. As a result, pollen analyses may underestimate the true rates of migration (i.e. as determined by first arrival followed by population expansion) in the case that a taxon is regionally rare though present at a site”.

As a result of these new analyses, we modified Fig. 2, Fig. 4, Fig. 5, Fig. 7, and the RMSE values in Table B3 (old Table B2), along with other parts in the manuscript text (refer to Comments 12, 14 and 20).

Comment 8: Table 1, Dispersal syndrome. Given that it is well established from observational data that Jay (Garrulus glandarius) collects, caches and then 'loses' significant numbers of beech nuts and acorns, I find it very strange (and inaccurate) that bird dispersal is not listed for Fagus and the various Quercus spp., for which taxa it is widely viewed as the most important mechanism of longer-distance dispersal. I am equally puzzled by the indication that Betula spp. and Picea abies are bird dispersed.  The only birds that I am aware of that feed upon the fruits/seeds of these taxa (e.g. Siskin - Carduelis spinus (Betula & Picea), Redpoll - Carduelis flammea/cabaret/hornemanni (Betula), Common Crossbill - Loxia curvirostra (Picea)) consume the fruits/seeds, usually on the tree, and do not, as far as I am aware, cache them or carry them away from the tree.  They may drop some fruits/seeds when feeding, but those will fall principally below the canopy of the tree.  The only basis that I can see for suggesting bird dispersal of these taxa is as an accidental mechanism, with the small winged fruits/seeds potentially adhering to damp feathers and thus being carried some distance from their source before being removed when the bird preens (Wilkinson, 1997).  I would consider such accidental dispersal as quite separate from a bird dispersal syndrome such as we see with Fagus and Quercus. It also seems to me that there is good evidence of dispersal of many nut-like fruits/seeds by larger mammals, the fruits/seeds being consumed but a proportion of them passing through the gut without digestion and being defecated (Giannakos, 1997).  Thus LA ought to be added for Corylus and Quercus spp.

Wilkinson, D.M. (1997) Plant colonization: are wind dispersed seeds really dispersed by birds at larger spatial and temporal scales? Journal of Biogeography, 24, 61-65.

Giannakos, P. (1997) Frugivory and seed dispersal by carnivores in the Rhodopi mountains of northern Greece. PhD, Durham.  (Available as an e-thesis to download.)

Response: Thank you for providing detailed information on the dispersal syndromes of our focal species. To be consistent, we run a second search on the TRY Database and applied occurrence threshold to define primary, secondary (<30%) and irrelevant (<5%) dispersal mechanisms. We further supported this search with information from research articles, including the articles you suggested in the comment. You may find details on species-specific dispersal syndromes and their references in the new Table B1. Accordingly, we updated Table 1, Dispersal syndrome.

Comment 9: L230: “LPJ-GM reduces the number of replicate units to one while using multiple explicitly placed patches per grid cell (1 km² each) in order to give a spatially explicit representation of the migration processes.”. Are the patches being simulated by LPJ-GM of 1 km²?  ...or is the half-degree grid cell being sub-divided into ~3000 (at the Equator) 1 km² grid cells, and then one LPJ-GM 'patch' (usually of extent 0.1 ha) being simulated for each 1 km² grid cell?  This requires clarification because, as written, the implication is that the 'patch' being simulated is of extent 1 km², which would be extremely computationally demanding, and hence is viewed as extremely unlikely.

Response: We clarified this part of the text (by also taking into account Major Comment 2 by the Anonymous Referee #1) as follows (L286-289): “The climate grid of 0.5 degrees resolution was subdivided into smaller cells of 1 km² area, where vegetation dynamics are simulated at a patch level of 1000 m², as is usually done for LPJ-GUESS simulations. Simulations were performed for a total of 500 years, covering an area of 201 x 201 cells with corridors located on the perimeter and the two major diagonals of the domain for a total of 1,197 simulated cells (see Fig. S1).”.

Comment 10: L233. As I noted in the Abstract, it would be preferable, and more accurate, to use 'propagule' rather than seed, throughout your text, as many of the dispersed entities are fruits (e.g. in the case of Betula).

Response: Here we decided to keep ‘seed’ as this term is used in reference to the model more than in an ecological sense. Please see also response to Comment 1.

Comment 11: L247-250: “This choice was based on evidences from pollen records (Birks and Birks, 2008) and phylogeographic studies (Palmé et al., 2003) suggesting a scenario of early colonization of treeless ground by Betula species within a very short time period (only hundreds of years) after the retreat of the ice sheet from northern Europe, followed by successive waves of colonization by later-successional tree species (Giesecke and Brewer, 2018).”. In many parts of northern Europe the first woody species to invade deglaciated areas were juniper (generally assumed to be Juniperus communis) and willow (usually assumed to be shrubby species of Salix). There is also the complication that where the evidence is from pollen analysis, Betula pollen could be derived from B. nana at least in the earliest stages of post-deglaciation colonisation. That said, by the time the first of the other tree species were expanding into these areas the woody vegetation was over large areas Betula woodland.  However, this was in many regions rapidly replaced by Pinus, presumably P. sylvestris, woodlands, so that later arriving temperate tree taxa were invading pinewoods, resulting in such no-analogue assemblages as pine-hazel-birch woodlands over moderately extensive areas in the early Holocene, that were subsequently invaded by temperate taxa.  Furthermore, I am not aware of any extensive areas where Fagus was migrating into Betula woodlands, not least because it reached northern Europe long after the birch woodland phase had passed; the same is true of Tilia cordata.  The westward expansion of Picea abies across Fennoscandia during the Holocene was largely into Pinus-dominated woodlands.... All in all, expansion into birch woodlands applies only to the minority of tree taxa that reached northern Europe during the first millennium of the Holocene.

Response: Thank you for mentioning other possible interactions with species other than Betula spp. in the early stage of colonization. Considering this information, we decided to include additional competitors (beside Betula spp.) for the relevant species. The text was modified as follows (L271-285): “This choice was based on evidences from pollen records (Birks and Birks, 2008) and phylogeographic studies (Palmé et al., 2003) suggesting a scenario of early colonization of treeless ground by Betula species, along with juniper (most likely Juniperus communis) and willow (as shrubby Salix species) (B. Huntley, 2022, pers. comm., January), within a very short time period (only hundreds of years) after the retreat of the ice sheet from northern Europe, and followed by successive waves of colonization by later-successional tree species (Giesecke and Brewer, 2018). Thus, the first tree species to expand after Betula spp. were likely invading large areas of birch woodland, though some regions were rapidly replaced by Pinus spp. (most likely P. sylvestris). For example, the westward spread of Picea abies across Fennoscandia was largely into Pinus-dominated woodlands. Furthermore, late-successional species such as Fagus and Tilia reached northern Europe after the phase of Betula-dominated woodlands had passed, and likely encountered other competitors in some parts of their expanding ranges (B. Huntley, 2022, pers. comm., January). Therefore, we decided to simulate the spread of Fagus sylvatica, Picea abies, and Tilia cordata with two additional competitors, i.e. Quercus robur, Pinus sylvestris, and Quercus robur or Pinus sylvestris, respectively (B. Huntley, 2022, pers. comm., January). Simulations with alternative competitors were conducted with optimized parameters (where the initial optimization was conducted assuming B. pendula as competitor) and evaluated in Sect. 3.3. This allowed assessing the impact of competition on the simulated migration rates for the three species.”. See also Comment 21.

Comment 12: L277: “Since we assumed maximum values of paleo-records obtained by classic pollen estimates to be generally over-estimations of true migration rates (see Sect. 2.1), we used the 75th percentiles of the observed range of migration values as the target observational value to improve and assess the performance of the model across species (see also Sect. 2.6.2).” See earlier comment; it is highly unlikely that the migration rates estimated from pollen data are over-estimates of true rates, for reasons that I have explained in that earlier comment.

Response: See response to Comment 7. We modified the text as follows (L325-328): “Since we assumed maximum values of paleo-records obtained by classic pollen estimates to be good estimates of the maximum potential spread achieved by species under optimal environmental conditions (and potentially under-estimations of true migration rates; see Sect. 2.1), we decided to use the upper boundaries of estimates from the literature (see Table 1) as the target observational value […]”.

Comment 13: Table 2, Parameter. All four of these parameters are going to have a positive influence on migration rate, albeit that the strength of their influence will vary enormously, with the dispersal distances having a vastly greater impact than fecundity or germination rate.  Indeed, if this is not what you find, then my conclusion would be that there was something wrong with your model!

Response: We agree. Our first way to validate the model was to check whether the relationship between migration rate and the four parameters followed the expected trend. We mentioned this in L263-266 (old L239-242).

Comment 14: L383-384: “We used the species-specific 75^th percentile of the observed values since upper boundaries (i.e. maximum values) are assumed to be potential over-estimations of post-glacial migration speeds (see Sect. 2.1).” That assumption again!  See earlier comments.

Response: See response to Comment 7. We modified the text as follows (L440-442): “We used the species-specific upper boundaries of the observed values since we assumed that the maximum rates detected by classic pollen analysis are good estimates of the maximum potential spread of a species under optimal climatic conditions (see Sect. 2.1).”.

Comment 15: L428-429: “while all parameters are overall positively related to migration rate over the range of each individual parameter”. This is no surprise; indeed, I predicted that this was the case when I saw the parameters listed in Table 2. It is also no surprise that migration rate is more sensitive to dispersal distances than to fecundity - this was something that Yvonne Collingham showed more than 25 years ago using her MIGRATE model (Collingham et al., 1996)

Collingham, Y.C., Hill, M.O. & Huntley, B. (1996) The migration of sessile organisms: a simulation model with measurable parameters. Journal of Vegetation Science, 7, 831-846.

Response: See response to Comment 13.

Comment 16: L447: “Picea abies (above 10 m yr^-1 from the average observed migration speed; Fig. 2…)”. Using a threshold of 10 m yr^-1 seems to me to be much too stringent, given the observed migration rate values in Table 1, and given the very considerable uncertainties in parameter estimates.

Response: Following your advice, we modified the threshold to exclude bad estimates from an absolute value (10 m yr^-1) to a percentage value of 15 % with respect to the maximum estimate of migration rate (i.e. simulated values 15 % below the maximum are under-estimations and vice-versa for over-estimations; 10 % was tested with identical results). Accordingly, the text was modified as follows: “Picea abies (<15 % from the maximum observed migration speed; Fig. 2…)” (L503-504).

Comment 17: L476-478: “the exponential power and twoDt functions generated good estimates for Fraxinus excelsior (RMSE = 25.75 % and 25.83 %, respectively), while the log-hyperbolic secant and logistic functions tended to over-estimate migration rates by 30-45 % for both species.” So far in this section Fraxinus excelsior is the only species that has been mentioned.  What is the other species being included here when you refer to 'both species'?

Response: Thank you for spotting this error. “for both species” was eliminated (L535).

Comment 18: L556: “when surpassing a certain biomass threshold (LAI = 0.5)…”. Surely LAI is not a 'biomass threshold'?  If you want to use a biomass threshold, then surely you should use the carbon-mass variable, rather than LAI, output by LPJ-GUESS?

Response: We changed “biomass threshold” to “leaf area index threshold” (L625). We used a LAI threshold as explained by Lehsten et al. (2019), https://doi.org/10.5194/gmd-12-893-2019.

Comment 19: L590-593: “In this case, the use of fat-tail kernels will produce a noisy and accelerating vegetation spread relative to a step-wise and slower spreading front given by exponentially-bounded kernels (Clark, 1998). This seems to suggest that occasional LDD events are more important for tree migration than local dispersal driven by more common vectors, at least in order to achieve the spreading rates of the paleo-records used in this study.”. Not only has the fundamental importance of rare long-distance dispersal events long been recognised, but we have demonstrated, using palaeoecological studies, that migration of at least one tree species during the Holocene progressed by long-distance dispersal events with subsequent infilling, rather than by a continuous 'migrating front'. Pollen analytical data will only rarely be able to show this, as Davis et al. (1991) demonstrated.  However, we were able to show this for the mid-Holocene expansion of Pinus sylvestris into the far north of Scotland by dendrochronological analyses of sub-fossil stumps buried in peat (Gear & Huntley, 1991). These recorded an initial colonisation event by very few individuals, followed, after a delay, by a second wave of establishment, presumably of progeny of the initial colonists, the delay probably reflecting the time taken for the initial colonists to reach an age sufficient to produce seeds.  This pattern was repeated at sites across the region, but with the initial colonisation taking place progressively later at sites more distant from the area occupied by P. sylvestris woodlands prior to the period of expansion (Huntley et al., 1997; Daniell, 1992, 1997).

Davis, M.B., Schwartz, M.W. & Woods, K. (1991) Detecting a species limit from pollen in sediments. Journal of Biogeography, 18, 653-668.

Gear, A.J. & Huntley, B. (1991) Rapid changes in the range limits of Scots Pine 4000 years ago. Science, 251, 544-547.

Huntley, B., Daniell, J.R.G. & Allen, J.R.M. (1997) Scottish vegetation history: The Highlands. Botanical Journal of Scotland, 49, 163-175.

Daniell, J.R.G. (1992) A dendrochronological study of subfossil Pinus sylvestris L. stumps from peat deposits at Badanloch in the far north of Scotland. M.Sc., University of Durham, Durham.

Daniell, J.R.G. (1997) The late-Holocene palaeoecology of Scots Pine (Pinus sylvestris L.) in north-west Scotland. Ph.D., University of Durham, Durham.

Response: Thank you for providing additional references to the topic. We added the following sentence in the Discussion (L662-666): “Supporting this, dendrochronological analyses of subfossil stumps have shown that the mid-Holocene expansion of P. sylvestris into northern Scotland likely progressed by LDD events with initial immigration by very few individuals followed by infilling after a delay (possibly by the mature progeny of the first colonizers), rather than by a step-wise migration front (Gear and Huntley, 1991; Daniell, 1992; Daniell, 1997; Huntley et al., 1997).”.

Comment 20: L634-635: “Pinus spp. have both the highest error relative to simulated migration rates and the largest uncertainty of observed migration rate across all species (from 600 to 1500 m yr^-1)”. Does this range of cited values really represent 'uncertainty'? - or do the different values reflect, as I strongly suspect, estimates made in different parts of the species' range boundary and/or during different stages of its range boundary expansion? Surely, the evidence from the estimate of 1500 m yr^-1 is that the species was able to achieve such a rate for at least some parts of its expanding range margin and at least at some time during its range expansion.  Thus, the target maximum rate for the model ought to be 1500 m yr^-1, and not any lower value (and not just 75% of this value!).

Response: See response to Comment 7. We removed the sentence “Pinus spp. have both…”.

Comment 21: L636-638: “We might have assumed the wrong competitor (see Sect. 2.3), e.g. some late-successional species (e.g. Fagus sylvatica) might have competed with additional species beside the early-successional Betula spp. (Giesecke and Brewer, 2018).”. See my earlier comment - you certainly did assume the incorrect competitor for many of the taxa investigated.  In the case of Fagus, a much more appropriate competitor across most of its range would be Quercus, which is a very different tree in many ways from Betula.

Response: We removed the above sentence as we have tested for alternative competitors as suggested (see Comment 11). Results from error and uncertainty analyses tested with alternative competitors (with respect to the default Betula spp.) showed overall small differences (L545-551): “Additional simulations with optimized parameters that tested alternative competitors (apart from Betula spp.) showed that migration rate is relatively robust to different competitors. For the three species tested with alternative competitors, the error between simulated migration speed and observed maximum value had a small increase, while staying below the threshold of good estimates (RMSE < 15 %; Sect. 2.5.2): Fagus sylvatica (from 5.3 % to 14 % with Quercus robur), Picea abies (from 1.8 to 11.8 % with Pinus sylvestris) and Tilia cordata (from 3.2 % to 3.6 % and 10.4 % with Pinus sylvestris and Quercus robur, respectively). Similarly, the overall error of the model had a minor increase (from 7.7 to 8.13 %) and the model utility stayed constant (0.58) (see Sect. 3.4).”.

Comment 22: L639-640: “We decided to simulate a homogeneous terrain (for permeability) and climate to ensure ideal climatic/topographic conditions for vegetation spread to match maximum migration speeds.” Arguably, it is the converse that it is true! Given that it is long-distance dispersal events that fundamentally drive the overall migration, and especially maintain the rapid rates, and also that in a real landscape such long-distance dispersal also is essential for most taxa if they are to disperse to new patches of suitable habitat (real landscapes never being homogeneous), I would argue that migration through a heterogeneous landscape is likely to be more rapid, as this is more likely to offer patches of suitable habitat that will support the species' migration.  Only in the (rare) case of a taxon for which a homogeneous area (which is in any case unrealistic) is uniformly suitable will such a landscape support rapid migration - but only for that taxon best suited to the uniform habitat offered.

Response: Thank you for your clarification. Our intent was to evaluate how our model simulated migration speed under “ideal” conditions (i.e. completely homogeneous terrain and non-limiting climate) so to reduce as much as possible the influence of climate and topography on the migration process, and concentrate on the effect of the model parameters. To clarify this point, we modified the sentence as follows (L706-711): “We decided to simulate the spread of tree species across a homogeneous terrain (for permeability) and non-limiting climate in order to reduce as much as possible the effect of the environment on the simulated migration rates. This allows assessing the impact of parameters and model structure, to conduct a comparison among species with equal simulation settings, and to optimize migration parameters to match maximum migration speeds from the literature (as we assume that a non-limiting climate may allow species to achieve the maximum potential spread recorder by paleo-records; see Sect. 2.1).”.

Comment 23: L641-642: “i.e. trees were mainly limited by dispersal and not by climate during the postglacial forest expansion of Europe.” I would argue that this issue was 'put to bed' by palaeoecologists more than 25 years ago. It was only 'resuscitated' following modelling by Svenning and colleagues (Svenning & Skov, 2004) using an unrealistic modelling approach that inevitably generated the results that they, sadly, were able to publish - an example where peer review was seriously ​deficient. Have a look at Prentice et al. (1991).

Svenning, J.C. & Skov, F. (2004) Limited filling of the potential range in European tree species. Ecology Letters, 7, 565-573.

Prentice, I.C., Bartlein, P.J. & Webb, T., III (1991) Vegetation and climate change in Eastern North America since the last glacial maximum. Ecology, 72, 2038-56.

Response: We eliminated this sentence.

Comment 24: L646: “…dispersal ability (e.g. adaptation to the new environment may reduce the need to disperse).” Has such adaptation been observed in (any) plants?  It seems to me to be extremely unlikely to develop in long-lived perennials such as trees over any timescale that humans might observe over. I am only aware of such adaptation being demonstrated in the case of some short-lived insects.  Whilst there is some evidence that populations of some bird species that undergo longer migrations have greater wing lengths than populations of the same species that are either shorter-distance migrants or residents, this is apparently a relatively long-term (at least multi-millennial) evolutionary adaptation rather than a short-term (decadal) adaptation to new circumstances, as far as I am aware.

Response: We clarified this point by adding a reference on dispersal adaptation for insects, whereas adaptation for trees is rare due to the slow response of long-lived perennials, such as trees (725-730): “Though there are some instances of evolutionary responses concerning dispersal ability (e.g. adaptation to the new environment may reduce the need to disperse), these are more prevalent among short-lived organism such as insects (De Meester et al., 2018). As long-lived perennials, trees seem to have a slow plastic/evolutionary response with respect to other taxa and are more likely to respond to climate change with range shifts (Berg et al., 2010; Lenoir et al., 2020; see also Sect 6.3 and Table 11 of Birks, 2019 for the little evidence of plasticity/adaptation of terrestrial plants during the Quaternary).”.

Comment 25: L662: “mid-Holocene estimates”. Not all the estimates are 'mid-Holocene', indeed a substantial proportion of them are from the earliest Holocene, with only the later migrating taxa, e.g. Fagus, Tilia, Picea, having estimates based largely on mid-Holocene migrations.

Response: We modified “mid-Holocene” to “early- and mid-Holocene”.

Comment 26: L664: “past and current global warming have different intensities and species are submitted to different conditions (e.g. more fragmented habitats and unlimited human-driven dispersal in present/future conditions”. Most habitats are naturally 'fragmented', in that they exist as discrete patches in a heterogeneous landscape.  What is much more important is that today many of those patches have been destroyed - i.e. it is habitat loss rather than fragmentation that is key. Furthermore, as Huntley & Collingham (2000) showed, a landscape with a given proportion of suitable habitat will better sustain migration rates if that fraction comprises a large number of small patches, as opposed to a smaller number of large patches.

Collingham, Y.C. & Huntley, B. (2000) Impacts of habitat fragmentation and patch size upon migration rates. Ecological Applications, 10, 131-144. Aside from quibbling about the idea that human-driven dispersal is 'unlimited', this has to be offset by the severe decrease in abundance, or even regional extirpation, of 'natural' agents of long-distance dispersal, especially larger mammals, but also birds, as a result of human activities, especially land-use and 'pollution' (where the latter includes pollution with biocides).

Response: We modified the text accordingly: “…to different conditions (e.g. habitat loss, decrease of natural dispersal vectors such as large mammals or birds, and increase of human-driven dispersal in present/future conditions with respect to the early- and mid-Holocene;” (L745-746) and “…trees will be limited by their dispersal ability (i.e. migration lag), and especially by LDD events in landscapes where the availability of suitable habitats has been severely decrease by human activities.” (751-753). Furthermore, we added a reference to the impact of habitat loss on species range shifts in the Introduction (L38-40): “Concerning the human impact on landscape reconfiguration, habitat loss has been shown to have a significant and negative impact on vegetation expansion both in models (Collingham and Huntley, 2000; Dullinger et al., 2015; Saltré et al., 2015) and real case studies in the current century (Guo et al., 2018).” (see also response to Comment 2 of the Anonymous Referee #2).

Comment 27: L667: “…we identified important mechanisms controlling migration rate (LDD events)”. As I have observed earlier, this is hardly a new finding.

Response: We considered the fact that previous models found the same relevant mechanism (LDD events) as a validation of our model structure (as reported in the Discussion). Here our main goal was to implement the most important mechanisms of the migration process according to our model in an efficient way. To clarify this point, we modified the sentence to “…we implemented important mechanisms controlling migration rate (LDD events) into the model structure (fat-tail kernels).” (L748-749) and in the Conclusions to “…the efficient inclusion of influential migration mechanisms (LDD events) in the model structure (via fat-tail kernels) can improve our confidence in range shift predictions…” (L764-765).

Comment 28: L670: “The high velocity of isothermal shift predicted for the 21st century suggests that trees will be limited by their dispersal ability (i.e. migration lag), and especially by LDD events in a fragmented landscape with few available establishment sites.” NO - in landscapes where the availability of suitable habitat patches has been very severely reduced by human activities.

Response: See response to Comment 26.

We would like to thank again all reviewers for taking the time to review our manuscript.