Articles | Volume 21, issue 5
Research article
18 Mar 2024
Research article |  | 18 Mar 2024

Frost matters: incorporating late-spring frost into a dynamic vegetation model regulates regional productivity dynamics in European beech forests

Benjamin F. Meyer, Allan Buras, Konstantin Gregor, Lucia S. Layritz, Adriana Principe, Jürgen Kreyling, Anja Rammig, and Christian S. Zang

Late-spring frost (LSF) is a critical factor influencing the functioning of temperate forest ecosystems. Frost damage in the form of canopy defoliation impedes the ability of trees to effectively photosynthesize, thereby reducing tree productivity. In recent decades, LSF frequency has increased across Europe, likely intensified by the effects of climate change. With increasing warming, many deciduous tree species have shifted towards earlier budburst and leaf development. The earlier start of the growing season not only facilitates forest productivity but also lengthens the period during which trees are most susceptible to LSF. Moreover, recent forest transformation efforts in Europe intended to increase forest resilience to climate change have focused on increasing the share of deciduous species in forests. To assess the ability of forests to remain productive under climate change, dynamic vegetation models (DVMs) have proven to be useful tools. Currently, however, most state-of-the-art DVMs do not model processes related to LSF and the associated impacts. Here, we present a novel LSF module for integration with the dynamic vegetation model Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS). This new model implementation, termed LPJ-GUESS-FROST, provides the ability to directly attribute simulated impacts on forest productivity dynamics to LSF. We use the example of European beech, one of the dominant deciduous species in central Europe, to demonstrate the functioning of our novel LSF module. Using a network of tree-ring observations from past frost events, we show that LPJ-GUESS-FROST can reproduce productivity reductions caused by LSF. Further, to exemplify the effects of including LSF dynamics in DVMs, we run LPJ-GUESS-FROST for a study region in southern Germany for which high-resolution climate observations are available. Here, we show that modeled LSF plays a substantial role in regulating regional net primary production (NPP) and biomass dynamics, emphasizing the need for LSF to be more widely accounted for in DVMs.

1 Introduction

In temperate climates, late-spring frost (LSF) plays a critical role in the functioning of forest ecosystems (Grossman2023). Below-freezing temperatures during the early stages of leaf development damage photosynthetic tissue (Inouye2000; Chen et al.2023), hamper secondary growth (Dittmar et al.2006; Príncipe et al.2017), induce mobilization of stored carbon reserves to repair damaged tissues (D'Andrea et al.2019), and ultimately constrain the range limits of affected tree species (Körner et al.2016; Kollas et al.2014a). Consequently, evidence suggests that LSF substantially reduces forest ecosystem productivity (Hufkens et al.2012).

The frequency of LSF has increased across Europe in recent decades, displaying a trend that will likely be exacerbated under a changing climate as an earlier start to the growing season leads the timing of leaf-out and periods with high likelihood of frost to increasingly coincide (Zohner et al.2020; Ma et al.2019). Due to increasing warming, many deciduous tree species have experienced a shift towards earlier budburst, flowering, and subsequent leaf development (Morin et al.2009; Menzel et al.2001). While a longer growing season facilitates forest productivity (Keenan et al.2014; Duveneck and Thompson2017), it may also bear a greater risk of detrimental LSF impacts as the period in which trees are most susceptible increases (Chamberlain and Wolkovich2021; Ma et al.2019; Sangüesa-Barreda et al.2021). Thus, it is unclear to what extent climate change will alter the impact of LSF on forest ecosystem productivity.

One of the main tree species in central Europe is Fagus sylvatica (European beech), which under potential natural vegetation conditions is the dominant species across large regions of the European landscape (Bohn and Welß2003). Although historic land management has reduced the proportion of European beech in forests, recent management efforts have focused on re-establishing a higher share of deciduous broadleaf species, including beech (Kenk and Guehne2001; Schütz1999). These efforts aim to increase forest resilience to climate change (Yousefpour et al.2018). In a simple twist of fate, European beech tends to be relatively susceptible to late-spring frost, which consequently co-determines its range limits across Europe (Gazol et al.2019; Menzel et al.2015; Kollas et al.2014a; Bolte et al.2009). Since beech has also been shown to be negatively impacted by drought events (Meyer et al.2020; Dulamsuren et al.2017; Zimmermann et al.2015; Scharnweber et al.2011), the success of recent forest transformation efforts hinges on the ability of beech forests to remain productive in their current distribution range under climate change.

In this context, dynamic vegetation models (DVMs) are a useful tool to assess the impact of climate change and extreme events on forest productivity (Yao et al.2022; Medlyn et al.2011; Gampe et al.2021; Rammig et al.2010). Nevertheless, most current state-of-the-art DVMs do not simulate LSF and the associated consequences, possibly overestimating the positive effects of climate change (e.g., longer growing seasons) on the carbon sequestration potential of temperate forest ecosystems (Liu et al.2018) and underestimating detrimental effects from climate-change-related impacts on tree health and mortality. Considering forests currently account for nearly 50 % of terrestrial net primary productivity (NPP) (Bonan2008), the ability of vegetation modeling to assess the impact of current and future LSF on forest NPP is crucial. To bridge this gap, we present a novel late-frost module for integration with the dynamic vegetation model LPJ-GUESS (Smith et al.2014; Hickler et al.2012). A phenological model to predict budburst introduced by Kramer et al. (2017) forms the basis of this new model and allows us to connect the start of leaf development with the occurrence of sub-zero temperatures. This integration with LPJ-GUESS – named LPJ-GUESS-FROST – enables us to model productivity dynamics and directly attribute them to the presence (or absence) of LSF.

The new model implementation is validated against a regional tree-ring network, and subsequently, we present the effects of LSF on European beech NPP and carbon stocks, determine shifts in intra-tree carbon allocation patterns as predicted by the model, and identify regional LSF hotspots across our study region.

2 Methods

In 1953 and 2011, well-documented late-frost events occurred across large regions of Bavaria. The 1953 event was centered around the Alps and alpine foothills in southern Bavaria; the 2011 event had its epicenter in Franconia in northern Bavaria. In both cases, freezing damage was observed in European beech (Príncipe et al.2017; Dittmar et al.2006). We make use of these two observed events to validate our frost implementation in LPJ-GUESS-FROST and analyze the impact of LSF on European beech productivity and biomass across Bavaria.

2.1 Overview of the dynamic vegetation model LPJ-GUESS

We use LPJ-GUESS (version 4.0.1; Lindeskog et al.2017), a well-established dynamic vegetation model designed to simulate ecosystem processes on a regional to global scale (Smith et al.2001, 2014). Vegetation is represented by plant functional types (PFTs), which cycle through establishment, growth, competition, and mortality. Generally, a PFT groups attributes (phenology, life strategy, drought tolerance, bioclimatic limits, etc.) of multiple similar individual species and represents them through a set of parameters. In this study, we follow the commonly used parameterization of Hickler et al. (2012) developed specifically for European tree species to explicitly simulate European beech. The model is driven by gridded daily climate, soil texture, nitrogen deposition, and global atmospheric CO2 inputs (for a more detailed description, see Sect. 2.6). Processes are modeled on a grid cell basis where the spatial resolution of the grid cells follows the spatial resolution of the climate inputs.

The processes modeled by LPJ-GUESS to simulate primary production and growth include photosynthesis and stomatal conductance based on BIOME3 (Sykes and Prentice1996); allocation of NPP to various compartments based on allometric constraints (Sitch et al.2003); stochastically simulated population dynamics (Hickler et al.2004); biomass-destroying disturbance; and nitrogen, soil, and litter processes (Smith et al.2014). Simulated vegetation dynamics emerge from the interaction between growth and competition for resources (e.g., light, water, and nutrients). For each grid cell, multiple “replicate” patches are simulated, each of which represents a random sample of the grid cell to account for idiosyncracies in disturbance and stand development of different vegetation stands. Within each patch, a single average individual represents cohorts of individuals established in a given year, with all individuals of a given cohort sharing the same size and form as they grow.

A key process in the context of this study is summergreen-leaf phenology, which is explicitly modeled in LPJ-GUESS as a function of daily mean temperature following a generalized chilling and growing-degree-day model (Smith et al.2001). Chill days occur when daily mean temperature falls below 5 °C and, combined with PFT-specific parameters regulating the thermal requirement for budburst (GDDbase; Eq. 1), determine the length of dormancy.

(1) GDD base = k chill , a + k chill , b × e k chill , k × chilldays

Subsequently, when a sufficient growing-degree sum (GDD5; Eq. 2, mean temperature (T) over 5 °C) for a given PFT has accumulated,

(2) GDD 5 = max { 0 , T - 5 } ,

leaf unfolding commences, and the phenological status (phen, Eq. 3) is updated daily as a fraction of complete canopy cover, ranging between 0 and 1, until maximum canopy cover is reached

(3) phen = 0 , GDD 5 GDD base min { 1 , GDD 5 - GDD base GDD ramp } , GDD 5 > GDD base ,

where GDDramp is a PFT parameter.

2.2 Implementation of budburst model and late-spring frost in LPJ-GUESS

We extended the leaf phenology calculation in LPJ-GUESS by implementing a novel model simulating LSF, running in parallel to the standard LPJ-GUESS phenology module (see Fig. A2). Late-spring frost is primarily a disruption of the phenological cycle of leaf development caused by the overlap of leaf-out and a critical sub-zero leaf temperature. The standard implementation of LPJ-GUESS phenology merely calculates the fraction of complete canopy cover (see above); thus, we here used a more specific model for calculating budburst status from Kramer et al. (2017). This sequential two-stage chilling and forcing model relies on daily mean temperature to calculate the budburst status, returning 0 when budburst has not yet occurred and 1 when budburst has taken place. The first stage of the model is the chill state (Sc; Eq. 4), which determines the period of phenological rest, calculated as

(4) S c ( t ) = i = t 0 t R c ( T i ) ,

where t0 is 1 November (the day on which the phenological model resets), t is the current time step (i.e., day), and Ti is the mean air temperature at time step i. Rc(Ti) is the rate of chilling as a function of the mean air temperature at time i, formulated as

(5) R c ( T ) = 0 , T < T c , min T - T c , min T c , opt - T c , min , T c , min T T c , opt T - T c , max T c , opt - T c , max , T c , opt T T c , max 0 , T > T c , max ,

where Tc,min, Tc,opt, and Tc,max are the minimum, optimum, and maximum mean air temperature required to advance chilling.

Analogously, the forcing state (Sf; Eq. 6) is calculated as the running sum of the forcing rate (Rf; Eq. 7).


where t1 is the time at which the chill state reaches the critical value Sc,crit (a species-specific constant derived from empirical observations determining the end of the chilling period), Tf,min is the minimum mean air temperature required for forcing, and af and bf are species-specific fitted constants. Subsequently, budburst occurs (i.e., model state is 1) when the forcing state (Sf) attains a critical value Sf,crit, which, like the critical value for chilling (Sc,crit), is derived from empirical observations.

In a second step, we cross-reference the budburst status (0 or 1) with the daily minimum temperature (Tmin) to determine late-frost status. When budburst has already occurred (budburst status is 1) and Tmin crosses the temperature threshold Tfrost, late frost occurs. Accordingly, we calculate LSF as

(8) LSF ( leafout , T frost , T min ) = 1 , leafout = 1  and  T min < T frost 0 , else .

To approximate the generally localized effect of LSF, we model Tfrost stochastically for each individual patch by randomly drawing from a Gaussian distribution with mean Tfrostµ and standard deviation Tfrostσ (for details, see Sect. 2.7 and Fig. A1). A full list of the parameter values for Eqs. (1)–(5) is shown in Sect. 2.7 (Table 1).

Late-frost damage is modeled as the phenology status (where 0 indicates no canopy and 1 indicates a full canopy) being reset to zero (i.e., complete removal of existing leaves), followed by a leafless period after which the phenological state continues to advance until it reaches full canopy coverage. We implemented a constant leafless period based on observations from several studies (Menzel et al.2015; Rubio-Cuadrado et al.2021; Nolè et al.2018), which found that the time between frost and full development of the second cohort of leaves ranged between  40 and  80 d. Consequently, we implement a 40 d leafless period where trees have no leaves, followed by a 20 d re-greening period where phenology steadily increases to 1 (i.e., full leaf coverage).

Since we model late-frost damage stochastically (as a function of Tfrostµ and Tfrostσ) for each patch, this means that any given grid cell will have late-frost damage ranging between 0 % and 100 %. For example, if none of the patches experience frost, the grid cell has 0 % frost damage. If half of the patches experience frost, the grid cell has 50 % frost damage. Through this mechanism, we approximate a continuous function for frost damage.

The additional carbon necessary to rebuild the canopy after LSF is drawn from the carbon storage pool in LPJ-GUESS-FROST. Carbon allocation in LPJ-GUESS occurs at the end of the year, by distributing the accumulated NPP to the various biomass compartments according to a set of allometric constraints. If the accumulated NPP is not sufficient to allocate enough carbon to each compartment to satisfy these constraints, additional carbon can be “borrowed” from a carbon storage pool. This pool must be refilled when enough NPP is available. The carbon lost due to LSF is calculated as the total carbon allocated to leaves for a given year multiplied by the fraction of canopy coverage at the time of LSF. At the end of the year, during allocation, this fraction of carbon is deducted from the storage pool. A portion of this is immediately “repaid”, reducing the amount of carbon available for allocation to the structural compartments. Over subsequent years, more NPP is allocated to the storage pool until it is completely refilled.

2.3 Site-level tree-ring data

To validate the results of LPJ-GUESS-FROST, we utilized a tree-ring network consisting of previously published data from 21 sites that were affected by LSF events in 1953 (Dittmar et al.2006) and 2011 (Príncipe et al.2017), supplemented by data from 5 sites at the epicenter of the 2011 LSF that have been newly sampled for this study. These data allowed us to retrospectively analyze the effect of LSF on productivity in European beech. Radial growth – expressed as ring width – integrates multiple signals aside from age and climate and can be used as an indicator of variation in forest productivity (Xu et al.2017). In the aggregate, the ring width of a given year is composed of age-, size-, climate-, and disturbance-related trends and additional, often unexplained, variability (Cook1987). Nonetheless, ring-width data have proven a useful tool to investigate the effect of climate on tree growth (e.g., Jevšenak2019; Anderegg et al.2020; Zang et al.2014; Bhuyan et al.2017; Wilmking et al.2020). To isolate the climate signal, all tree-ring-width data were detrended using a cubic spline with a frequency cutoff of 0.5 at 30 years. Detrending removes age-related growth trends from tree-ring data. The residuals of the determined spline, also called ring-width indices (RWIs), consequently depict mainly climate-induced growth variations (Sullivan et al.2016; Cook and Peters1997). Following detrending, the individual tree-ring series (ranging from 10 to 30 per site) were aggregated to site-level chronologies. The median, minimum, and maximum lengths of raw tree-ring chronologies are 113 years, 32 years, and 227 years, respectively. Since the climate data used to drive LPJ-GUESS and LPJ-GUESS-FROST begin in 1951, any RWI before that year were not included in subsequent analyses. Consequently, the median, minimum, and maximum lengths of the detrended ring-width chronologies used for comparison with LPJ-GUESS-(FROST) are 62 years, 33 years, and 68 years, respectively. The median, minimum, and maximum sample sizes for the site-level chronologies are 20, 10, and 30 individual tree-ring series, respectively (Fig. A3).

2.4 Forest condition monitoring data

To complement the validation of LPJ-GUESS-FROST via tree rings, we conducted an analysis of the 2011 LSF using remote-sensing data from the European Forest Condition Monitor (EFCM; Buras et al.2021). The EFCM provides high-resolution relative measures of forest greenness based on the Terra MODIS normalized difference vegetation index (NDVI) (Buras et al.2021). The key value used in the context of this study is the EFCM's proportional deviation from the median (PDM) metric. This metric facilitates the comparison of absolute changes in NDVI between pixels; i.e., pixels with a lower PDM can be said to display less relative greenness than pixels with a high PDM, as referenced against the long-term median (Buras et al.2021). To quantify the impact of the LSF on 4 May 2011, we used the PDMs from 25 May. Since the EFCM provides data integrated at a 16 d interval, this ensured that the entire post-LSF response was depicted in the PDMs. PDMs from an earlier time step would have necessarily also included pre-LSF signals. Additionally, we cross-referenced the PDMs with the regions of Bavaria that experienced sufficiently negative temperatures for LSF (see Table 1) on 4 May.

2.5 Climate data

We use two separate climate datasets to force LPJ-GUESS-(FROST) for this study. To reproduce the known site-level frost events in 1953 and 2011, we used historic, thin-plate spline-interpolated daily climate station observations (mean temperature, minimum temperature, and precipitation) from the German Meteorological Service (DWD), as downloaded from the climate data center (CDC;, last access: 15 December 2023). For each day over the analyzed period (1951–2020), available observations from an average of 227 climate stations (range 182–243 stations per day) were mapped to a digital elevation model (SRTM 90;, last access: 15 December 2023) with a spatial resolution of 250 m × 250 m using a 3-dimensional thin-plate spline (function TPS in the “fields” package) with longitude, latitude, and elevation as predictors for the temperature field. This was done to account for elevational effects on minimum temperatures, which are particularly important aspects when applying critical temperature thresholds for LSF mapping. Grid sizes of commonly used gridded products (e.g., 0.1°) render gridded products too coarsely in regions with a heterogeneous topography to represent small-scale variations in minimum temperature. The mean RMSE of mapped vs. observed temperatures was 0.43 ± 0.24 K (mean ± SD).

To analyze the regional variation in LSF dynamics across Bavaria, we instead used the BayObs gridded daily climate data product provided by the Bavarian Environment Agency (LfU). This dataset provides historical daily mean temperature, daily minimum temperature, and the daily precipitation sum from 1951 until 2020 on a 5 km spatial resolution for the Bavarian domain (Bayerisches Landesamt für Umwelt2020). We used this dataset to avoid the excessive computational load associated with running Bavaria-wide simulations at a resolution of 250 m by 250 m.

2.6 Modeling protocol

We conducted two separate simulation experiments. Firstly, to determine the ability of LPJ-GUESS-FROST to reproduce the effect of known LSF on European beech growth, we forced LPJ-GUESS and LPJ-GUESS-FROST with the historic spline-interpolated climate data for the 26 sites in our tree-ring network. Secondly, to ascertain the wider impact of and regional variation in LSF on European beech productivity and biomass dynamics, we forced LPJ-GUESS and LPJ-GUESS-FROST with the BayObs data for the entire Bavarian domain. Aside from the climate inputs, we used the same modeling protocol for both sets of simulation experiments.

To ensure that the simulated ecosystems were in near equilibrium at the start of the simulation experiment period, we spun up the model for 1000 years using recycled climate data from the first 30 years of the climate inputs. To establish a reference baseline, we used LPJ-GUESS for the spinup, that is, the version of the model without LSF. During the spinup, stochastic disturbances were turned on to facilitate a heterogenous age structure of the simulated forests. Subsequently, we used the identical post-spinup state to start the simulation experiment period (1951–2020) runs for both LPJ-GUESS and LPJ-GUESS-FROST. To isolate the effect of LSF on productivity and biomass dynamics, we switched off stochastic disturbances in the simulation experiment period. This ensured that any diverging responses in productivity were introduced only in the simulation period and allowed us to more accurately attribute the effect of LSF on these responses.

2.7 Model parameterization

We adapted the commonly used parameterization of LPJ-GUESS for the European tree species in Hickler et al. (2012) to simulate monospecific European beech stands. The central parameters for the new frost module are described in Table 1. We used 25 replicate patches and, when applicable (i.e., during spinup), a disturbance interval of 200 years.

The frost threshold, i.e., the temperature at which leaves of European beech are damaged, is a source of some uncertainty (Chen et al.2023). Commonly, a threshold of 2.2 °C is used for the false-spring index (Schwartz1993; Schwartz et al.2006; Ma et al.2019), and remote-sensing observations have indicated significant differences in canopy greenness between frost-affected and frost-unaffected beech at a minimum temperature of 1 °C (Buras et al.2021). Here, we conducted a sensitivity analysis to determine the frost threshold (1.65 °C ± 0.85 °C) by varying the frost threshold and assessing at which point the simulated response no longer matched the observed RWI response (see Fig. A1 for details). Our threshold falls well within the midrange of temperatures at which significant correlations of late-spring frost severity and gross primary productivity (GPP) were found by Chen et al. (2023).

Table 1Parameter values used for budburst and LSF determination in LPJ-GUESS-FROST. Values for budburst calculation (b) follow those published in Kramer et al. (2017). The parameter values for the frost threshold (a) were determined using a sensitivity analysis comparing LPJ-GUESS-FROST output with RWI data from our study sites (see Fig. A1 in the Supplement for details).

Download Print Version | Download XLSX

2.8 Model evaluation

To assess the efficacy of LPJ-GUESS-FROST at simulating LSF, we cross-referenced simulated NPP with the observed RWI. For both metrics (NPP and RWI), we computed the resistance index (Rt) introduced by Lloret et al. (2011). This index quantifies the ratio of radial growth in an event year (e.g., a frost year) to pre-disturbance growth. In this context, growth during and before a disturbance are defined as the average growth performance across a fixed time period (Lloret et al.2011). The direct impacts of LSF dynamics are confined to a single vegetation season; consequently, we only consider growth performance in the frost year to quantify performance during disturbance. We used a 2-year pre-disturbance window – in contrast to the more commonly used 3-year window (Pretzsch et al.2013) – since the first LSF occurred in 1953 and our climate data start in 1951. Hence, Rt in this study solely refers to the impact of LSF on growth in comparison to the 2 years preceding LSF. We assessed statistical differences in the Rt to LSF of LPJ-GUESS, LPJ-GUESS-FROST, and RWI using a pairwise Wilcoxon rank test (Bauer1972).

Additionally, to analyze the impact of LSF on European beech productivity across Bavaria, we calculated the loss of NPP and biomass due to LSF as the difference between the output of LPJ-GUESS and LPJ-GUESS-FROST during the post-spinup simulation period from 1951–2020.

Post-processing of model output, subsequent statistical analysis, and manuscript authoring were done in R version 4.2.1 (R Core Team2022) with the addition of the following packages: meta-package tidyverse (Wickham et al.2019), dplR (Bunn et al.2022), fields (Nychka et al.2021) ggthemes (Arnold2021), here (Müller2020), janitor (Firke2021), multcompView (Graves et al.2019), ncdf4 (Pierce2022), patchwork (Pedersen2020), RANN (Arya et al.2019), rnaturalearth (Massicotte and South2023), rnaturalearthhires (South et al.2024), terra (Hijmans2022), scico (Pedersen and Crameri2022), sf (Pebesma2018), and zoo (Zeileis and Grothendieck2005).

3 Results

3.1 Model validation

The negative impact of LSF on RWI is evident in both documented frost years. In 1953, RWI was reduced to nearly 45 % of the pre-frost period. The impact of the 2011 frost was less severe, reducing RWI to roughly 70 % of the pre-frost baseline (Fig. 1).

LPJ-GUESS simulated a slightly higher NPP in both frost years than in the pre-frost period (indicated by Rt>1 in Fig. 1), and in both 1953 and 2011 Rt from LPJ-GUESS was significantly different from RWI Rt. In contrast, LPJ-GUESS-FROST simulated a substantial reduction in NPP (Rt<1) in both years. In 2011, the impact of the LSF on productivity is of a similar magnitude in both the tree-ring data (Rt̃=0.69) and LPJ-GUESS-FROST (Rt̃=0.55). For the 1953 LSF, the pattern is similar. LPJ-GUESS-FROST simulates a significant NPP reduction (Rt̃=0.57) in response to LSF, and no significant difference is seen between LPJ-GUESS-FROST and the tree-ring data (Rt̃=0.45).

For the 2011 LSF, the remote-sensing analysis confirms the negative impact due to frost. In the areas of Bavaria that experienced temperatures below the frost threshold of 1.65 °C on 4 May, the PDMs were substantially lower than in the areas not affected by LSF (Fig. A4).

Figure 1Resistance of observed radial growth and simulated NPP to two well-documented late-frost events. LPJ-GUESS-FROST manages to capture the late-frost signal in both cases, albeit with a larger spread in 1953. The boxplots show the median, quartiles, and 1.5 interquartile range. Lettering indicates homogenous groups based on a pairwise Wilcoxon rank sum test with Bonferroni correction for multiple testings.


Figure 2Time series of annual net primary productivity (NPP) from LPJ-GUESS and LPJ-GUESS-FROST. The solid lines show the mean NPP across all 2866 grid cells. The shaded areas contain the 95th percentile range of NPP. Note that due to lower NPP caused by LSF in LPJ-GUESS-FROST biomass is lower, leading to a decrease in maintenance respiration. In non-frost years, the effect of lower maintenance respiration can lead to higher NPP in LPJ-GUESS-FROST than in LPJ-GUESS.


3.2 Effects on productivity

The range of NPP across all grid cells in Bavaria varied from nearly 0.3 kg C m−2 to around 0.6 kg C m−2 in LPJ-GUESS. Introducing late-frost dynamics increased the variation in NPP to range from ca. 0.15 to 0.6 kg C m−2 across all grid cells as some regions suffered from heavily decreased productivity in response to late-frost damage (Fig. 2). Simulated NPP in frost years was roughly 50 % lower than in non-frost years. Averaged across the entirety of Bavaria, the cumulative reductions in NPP caused by LSF resulted in 0.79 kg C m−2 of lost net productivity by the end of the simulation period in 2020 (Fig. 3a).

Figure 3(a) Cumulative differences in NPP between LPJ-GUESS-FROST and LPJ-GUESS due to late-frost impacts. The red line represents the mean cumulative NPP difference across all 2866 grid cells. The shaded area contains the 95th percentile range of cumulative NPP differences. (b) Differences in total vegetation carbon biomass due to late-frost impacts. The red line shows the difference in mean carbon mass (in kg m−2) between LPJ-GUESS-FROST and LPJ-GUESS across all 2866 grid cells. The shaded area contains the 95th percentile range of carbon mass. (c) Depiction of each tissue pool's contribution to the total difference in vegetation carbon biomass between LPJ-GUESS-FROST and LPJ-GUESS.


The lost productivity translates to biomass loss. By 2020, the effects of LSF resulted in a mean loss of 0.57 kg m−2 in vegetation carbon. For the 95th percentile of simulated grid cells, the change in vegetation carbon by the end of the simulation in 2020 ranges from a loss of 1.22 kg m−2 to a gain of 0.07 kg m−2 (Fig. 3b). This biomass loss primarily affects the sapwood, which accounts for 0.38 kg m−2 of lost vegetation carbon by 2020 (Fig. 3c).

The modeled biomass loss is in agreement with regional altitudinal patterns across Bavaria: regional maxima of biomass loss are concentrated in low-mountain areas in the south (alpine foothills), southeast (the Bavarian Forest), and northern parts (Franconian Jura) of Bavaria (Fig. 4).

Figure 4(a) Map showing the regional effects of late frost on whole-tree biomass loss by the end of the simulation period (2020). The heaviest losses, with up to  25 % loss of biomass due to LSF, were simulated in low-mountain regions across Bavaria (alpine foothills, the Bavarian Forest, and Franconian Jura) as indicated by (b) the elevation map (based on data from Copernicus DEM; European Space Agency2024). (c) The highest biomass losses (expressed in percentage of loss compared to the simulation without frost) tended to occur at higher elevations. On the contrary, low biomass losses are more evenly spread across all elevations.

4 Discussion

Using our novel implementation of late-spring frost dynamics in LPJ-GUESS (termed LPJ-GUESS-FROST), we managed to simulate the observed effect of two distinct frost events across large regions of Bavaria (Fig. 1). In both LSF years (1953 and 2011), the results from LPJ-GUESS-FROST match the RWI-based observations of variability in productivity in showing the distinct negative impact of LSF. There are no significant differences between the model and observations for both the 1953 and 2011 LSF. Notably, however, LPJ-GUESS-FROST simulates a more heterogenous frost response in 1953 than in 2011. This residual difference in the responses of LPJ-GUESS-FROST in 1953 and in 2011 can be explained by the simulated onset of budburst in those years. In 2011, the simulated onset of budburst across the 14 sites that experienced LSF ranged from day-of-year (DOY) 111 (21 April) to DOY 116 (26 April), well before LSF occurred between DOY 123 (3 May) and DOY 125 (5 May). Subsequently, from a phenological perspective, trees at all 14 sites were at risk of frost damage between DOY 123 and DOY 125. In contrast, in 1953 the simulated budburst across the 12 sites affected by LSF varied across a larger range, from DOY 117 (27 April) to DOY 142 (22 May). The recorded LSF took place between DOY 128 (8 May) and DOY 131 (11 May). At 3 sites, simulated budburst occurred after DOY 131, meaning that of the 12 sites, only 9 were phenologically predisposed to frost damage in 1953.

While the modeled responses to LSF are not significantly different from the observations, some differences remain. Discrepancies between the model output and the observational data (e.g., differences in range of modeled and observed responses in 2011) are to be expected. Firstly, we used gridded climate data to drive LPJ-GUESS-FROST, which almost certainly does not capture the actual local temperatures experienced by the sampled trees in either 1953 or 2011 as measured 2 m air temperature often deviates from canopy temperature during cold, clear nights (Kollas et al.2014b). Secondly, while ring-width indices have been shown to be a good proxy for annual variation in NPP, some mismatch between the two metrics must be expected (Xu et al.2017) due to differences in the demographic structure between modeled and observed tree stands. Nevertheless, our results demonstrate the efficacy of LPJ-GUESS-FROST at simulating the real-life impact of LSF on European beech productivity. Aside from reproducing the impact of the two well-known frost events for which tree-ring data were available, LPJ-GUESS-FROST simulates several additional LSF across Bavaria, suggesting that LSF is not a rare phenomenon in beech forests. Anecdotal evidence from agricultural observations confirms the occurrence of LSF in a number of the years simulated as LSF years by LPJ-GUESS-FROST (e.g., 2014, 2016, and 2019;, last access: 15 December 2023,, last access: 15 December 2023).

Consequently, LPJ-GUESS-FROST enables us to show the potential extent of losses in productivity and biomass due to frost damage (Fig. 3). Frost damage consistently led to lower productivity across large regions of Bavaria. Roughly 30 % of GPP and nearly 50 % of NPP were lost to frost damage in years with late frost. This matches results from Urbanski et al. (2007), who found a severe anomaly in net ecosystem exchange (40 % of the 13-year mean) following a late-frost event in Harvard Forest in 1998. Similarly, remote-sensing analysis of an LSF event in the northeastern USA in 2010 indicated a 7 %–14 % decrease in gross ecosystem productivity due to frost damage (Hufkens et al.2012). It is important to note that the remote-sensing approach used by Hufkens et al. (2012) necessarily includes information on all tree species in the study region. On the other hand, our simulations were specifically tailored to identify the impact of late frost on a single species, Fagus sylvatica. Therefore, a mismatch in productivity losses between the two studies is to be expected. The regional assessment done by Hufkens et al. (2012) included an early-leafing species, sugar maple, and two species with later leaf-out, American beech and yellow birch. Accordingly, sugar maples were most affected by frost damage. The lesser-affected species may in turn have buffered some of the productivity response to frost. These dynamics are absent from our study as we focused solely on simulating monospecific European beech stands. While this does not represent the current state of Bavarian forests, which until recently have been heavily managed to favor coniferous species for wood production (Kenk and Guehne2001; Schütz1999), historically, central European forests were dominated by beech (Ellenberg et al.2010). As management efforts increasingly aim to reinstitute large shares of beech (Kenk and Guehne2001), our aim is to highlight the potential effect of LSF on productivity dynamics in beech forests.

Additionally, we were able to simulate the extent to which losses in primary productivity translate to losses in tree biomass (Fig. 3). The majority of simulated biomass losses stem from reduction in sapwood biomass. This behavior is consistent with observed late-frost damage in tree rings where frost-damaged trees displayed smaller sapwood increments than their non-damaged counterparts (Rubio-Cuadrado et al.2021). Frost damage primarily manifests as a disruption of the photosynthetic apparatus via partial or full canopy defoliation (Menzel et al.2015; Inouye2000). Following defoliation, affected trees must recover their canopy before full photosynthetic activity can resume, effectively shortening their growing season (Augspurger2009). The additive effects of a shorter growing season and reallocation of stored reserves (D'Andrea et al.2019) to the new canopy consequently contribute to reduced radial growth in frost years. We capture part of this process by implementing a leafless period after LSF in LPJ-GUESS-FROST. After late-frost-induced canopy defoliation occurs in the model, the simulated phenology remains dormant for an extended period (see Menzel et al.2015; Rubio-Cuadrado et al.2021; Nolè et al.2018). The absence of leaves in the model prohibits photosynthesis and ultimately reduces the simulated annual NPP consistent with observations.

While we are able to compare our simulated productivity losses with those estimated by Hufkens et al. (2012) and Urbanski et al. (2007), to the best of our knowledge no previous study has quantified the effect of late frost on tree biomass. This would require forest stands with identical environmental conditions aside from late frost, which is not possible in a natural setting but can be simulated using DVMs. We show that the impact of LSF on carbon storage in plant biomass in European beech forests is non-negligible. In the final year of our 69-year simulation period, the difference in vegetation carbon between LPJ-GUESS and LPJ-GUESS-FROST amounted to a 5 % reduction due to the effects of LSF. To put this in context, Lindeskog et al. (2021) found that accounting for thinning in LPJ-GUESS yielded a reduction in vegetation carbon of 3 % to 5 % across Europe until 2010 compared to a simulation without thinning, suggesting that the effect of LSF on tree biomass can have a substantial effect on forest structure and functioning.

Nevertheless, these results must be interpreted with some caution. The implementation of late-spring frost in LPJ-GUESS-FROST is intended to provide potential estimates of frost-induced shifts in carbon dynamics. In reality, LSF is a highly localized disturbance, dependent on a forest stand's microclimate and frost tolerance. Accounting for these two factors introduces stochasticity into the frost scheme of LPJ-GUESS-FROST. Currently, LPJ-GUESS-FROST cannot simulate microclimate. While replicate patches are used to abstract forest structure, this is not the case for climate, which is constant across all patches in a grid cell. In addition, the frost tolerance of European beech leaves is a point of contention. Studies directly applying freezing temperatures to twig samples have found frost tolerance for European beech to range from 4.8 to 6.4 °C at and directly following budburst (Lenz et al.2013, 2016; Vitra et al.2017). These results conflict with observations of frost damage in European beech stands where recorded ambient temperatures from nearby climate stations were significantly higher. In fact, frost damage in European beech has been observed at temperatures as high as 1.2 °C (Príncipe et al.2017), and evidence from remote-sensing of LSF suggests canopy decline may already occur at 1 °C (Buras et al.2021). Similarly, 2.2 °C is a commonly accepted species-agnostic temperature threshold for palpable frost risk (Chamberlain et al.2019; Schwartz1993). This discrepancy is caused by the effect of radiative cooling. During clear windless nights, the temperature at the leaf tissue, where frost damage occurs, can potentially be several degrees lower than the ambient air temperature measured at climate stations (Matsui et al.1981; Neuner2014). To overcome these problems, we inverted the mechanism determining frost occurrence. Since we cannot meaningfully model microclimate on a patch level, and there is uncertainty in the specific frost threshold of leaves, we instead decided to model the frost threshold as stochastically variable at the patch level. In this manner we approximate the local differences in frost occurrence, while simultaneously accounting for the lack of an accurate frost threshold. While this abstraction approximates the real-world heterogeneity inherent to LSF, it is stochastic in nature. Therefore, the intensity of simulated frost damage may not always match observations.

5 Conclusions and outlook

Our results indicate the importance of including late-spring frost dynamics in DVMs. We demonstrate that known patterns of productivity loss in European beech due to LSF can be reproduced by LPJ-GUESS-FROST. Additionally, we found that LSF leads to distinct regional variation in simulated biomass loss and influences the allocation of carbon within individual trees. While these findings are relevant in their own right, they also imply a need to focus future research efforts on identifying the implications of LSF across multiple tree species and climate change scenarios. To do this, further modeling efforts in the realm of LSF should attempt to simulate microclimate within stands or, at the very least, develop routines to accurately account for the discrepancy between measured 2 m air temperature and leaf temperature due to radiative cooling. Considering that commonly used forcing from general circulation model (GCM) outputs is quite coarse and therefore unlikely to capture such localized effects, this poses a particular challenge for modeling future dynamics. Additionally, due to the strong controls of leaf-out on frost risk (Vitasse et al.2014; Ma et al.2019), improving the phenological models used in LPJ-GUESS-FROST should be a priority. Aside from parameterizing the phenological model for a wide range of broadleaf deciduous species (maple, hornbeam, ash, etc.), future phenology routines should also consider age-dependent variation in budburst times. Similarly, while we demonstrate that accounting for loss of photosynthetic capacity due to defoliation can simulate productivity losses that are comparable to observations, the relatively simple representation of carbon costs related to LSF damage and subsequent canopy rebuilding leaves room for improvement. Here, the integration of explicitly modeled non-structural carbohydrates into LPJ-GUESS-FROST could pave the way for a better representation of LSF impacts.

Appendix A

A1 Frost threshold sensitivity analysis

Figure A1Sensitivity analysis to determine the range of the frost threshold. For both known frost years, we ran simulations for LPJ-GUESS-FROST with varying frost threshold temperatures from 0.3 to 5.5 °C. We then computed the resistance index (see Methods) for each simulation and for the RWI-based observations. To determine Tfrostµ and Tfrostσ, we discarded all frost thresholds which resulted in a median resistance greater than 1 (i.e., higher productivity in a frost year than in pre-frost years) and computed the mean and standard deviation of the remaining thresholds (mean = 1.65; standard deviation = 0.85). The boxplots show the median resistance, quartiles, and 1.5 interquartile range. The solid red line indicates a resistance of 1.


A2 Schematic of phenology modules in LPJ-GUESS-FROST

Figure A2Schematic showing the integration of the new frost modules (red, green, and yellow) with the existing LPJ-GUESS phenology model (grey) and dependence on climatic drivers (blue). The LPJ-GUESS growing-degree day model (Eqs. 1–2 in the Methods section) simulates continuous leaf development from 0 (no leaves) to 1 (full canopy cover). This phenological status is then used for further model processes such as photosynthesis. The new frost module includes a parallel phenological budburst model which simulates a distinct point of budburst (i.e., the point at which developing leaves become sensitive to LSF). This model uses a sequential two-stage approach with a chilling stage (Rc, Sc) and a forcing stage (Rf, Sf) described in Eqs. (3)–(6) in the Methods section. This budburst status is used to determine LSF damage in conjunction with the minimum temperature (Tmin). In the case of LSF, the continuous phenological status is reset to 0, and a second cohort of leaves has to be rebuilt before photosynthetic activity can resume.


A3 Site-specific tree-ring chronologies

Figure A3Spline-detrended tree-ring chronologies per site. The solid black line depicts the mean chronology of ring-width indices (left y axis), and the grey shaded area signifies the number of individual tree-ring series contained in the mean chronology in any given year (right y axis). The dashed vertical lines indicate the LSF in 1953 and 2011. Superscripts indicate the reference for each site: 1 Meyer et al. (2020), 2 Meyer (2024c), 3 Dittmar et al. (2006), and 4 Principe et al. (2017).


A4 Geographical extent of negative PDMs indicating LSF impacts in 2011

Figure A4Map of proportional deviations from the median for the 2011 LSF from the EFCM. Black borders indicate the areas which experienced temperatures below 1.65 °C on 4 May 2011. The black dots show the tree-ring sites affected by LSF in 2011.

Code and data availability

LPJ-GUESS version 4.0.1 is openly accessible and available at (Lindeskog et al.2017). The modified version of LPJ-GUESS used in this study is openly accessible and available at (Meyer2024a). The data analysis pipeline to reproduce the results and figures from this study is openly accessible and available at (Meyer2024b). All data needed to reproduce the analysis are openly accessible and available at (Meyer2024c).

Author contributions

BFM, CSZ, and AB initiated the study. BFM performed the data analysis and model runs and prepared the manuscript. CSZ, AB, AP, and JK contributed data. KG, LSL, and AR contributed input to the model development. All authors contributed input to paper writing and revisions.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Biogeosciences. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We sincerely thank Christoph Dittmar for providing and allowing us to publish the tree-ring data from Dittmar et al. (2006) used in this study. We thank Anja Žmegač for preparing, dating, and measuring the newly sampled tree cores used in this study.

Financial support

This research has been supported by the Bayerisches Staatsministerium für Wissenschaft und Kunst (bayklif).

Review statement

This paper was edited by Lutz Merbold and reviewed by two anonymous referees.


Anderegg, W. R. L., Trugman, A. T., Badgley, G., Konings, A. G., and Shaw, J.: Divergent Forest Sensitivity to Repeated Extreme Droughts, Nat. Clim. Change, 10, 1091–1095,, 2020. a

Arnold, J. B.: ggthemes: Extra themes, scales and geoms for “ggplot2”, (last access: 3 March 2023), 2021. a

Arya, S., Mount, D., Kemp, S. E., and Jefferis, G.: RANN: Fast nearest neighbour search (wraps ANN library) using L2 metric, (last access: 3 March 2023), 2019. a

Augspurger, C. K.: Spring 2007 Warmth and Frost: Phenology, Damage and Refoliation in a Temperate Deciduous Forest, Funct. Ecol., 23, 1031–1039,, 2009. a

Bauer, D. F.: Constructing Confidence Sets Using Rank Statistics, J. Am. Stat. Assoc., 67, 687–690,, 1972. a

Bayerisches Landesamt für Umwelt [Hrsg.]: Bayerische Klimadaten – Beobachtungsdaten, Klima – Projektionsensemble Und Klimakennwerte Für Bayern, (last access: 7 May 2022), 2020. a

Bhuyan, U., Zang, C., and Menzel, A.: Different Responses of Multispecies Tree Ring Growth to Various Drought Indices across Europe, Dendrochronologia, 44, 1–8,, 2017. a

Bohn, U. and Welß, W.: Die Potenzielle Natürliche Vegetation, Klima, Pflanzen-Und Tierwelt, in: Leibnitz-Institut für Länderkunde, Nationalatlas Bundesrepublik Deutschland, 3, 84–87, 2003. a

Bolte, A., Ammer, C., Löf, M., Madsen, P., Nabuurs, G.-J., Schall, P., Spathelf, P., and Rock, J.: Adaptive Forest Management in Central Europe: Climate Change Impacts, Strategies and Integrative Concept, Scand. J. Forest Res., 24, 473–482,, 2009. a

Bonan, G. B.: Forests and Climate Change: Forcings, Feedbacks, and the Climate Benefits of Forests, Science, 320, 1444–1449,, 2008. a

Bunn, A., Korpela, M., Biondi, F., Campelo, F., Mérian, P., Qeadan, F., and Zang, C.: dplR: Dendrochronology Program Library in R,, (last access: 24 January 2024), 2022. a

Buras, A., Rammig, A., and Zang, C. S.: TheEuropean Forest Condition Monitor: Using Remotely Sensed Forest Greenness to Identify Hot Spots of Forest Decline, Front. Plant Sci., 12, 689220,, 2021. a, b, c, d, e

Chamberlain, C. J. and Wolkovich, E. M.: Late Spring Freezes Coupled with Warming Winters Alter Temperate Tree Phenology and Growth, New Phytol., 231, 987–995,, 2021. a

Chamberlain, C. J., Cook, B. I., García de Cortázar-Atauri, I., and Wolkovich, E. M.: Rethinking False Spring Risk, Glob. Change Biol., 25, 2209–2220,, 2019. a

Chen, L., Keski-Saari, S., Kontunen-Soppela, S., Zhu, X., Zhou, X., Hänninen, H., Pumpanen, J., Mola-Yudego, B., Wu, D., and Berninger, F.: Immediate and Carry-over Effects of Late-Spring Frost and Growing Season Drought on Forest Gross Primary Productivity Capacity in the Northern Hemisphere, Glob. Change Biol., 29, 3924–3940,, 2023. a, b, c

Cook, E. R.: The Decomposition of Tree-Ring Series for Environmental Studies, Tree-Ring Bull., 47, 37–59, 1987. a

Cook, E. R. and Peters, K.: Calculating Unbiased Tree-Ring Indices for the Study of Climatic and Environmental Change, Holocene, 7, 361–370,, 1997. a

D'Andrea, E., Rezaie, N., Battistelli, A., Gavrichkova, O., Kuhlmann, I., Matteucci, G., Moscatello, S., Proietti, S., Scartazza, A., Trumbore, S., and Muhr, J.: Winter's Bite: Beech Trees Survive Complete Defoliation Due to Spring Late-Frost Damage by Mobilizing Old c Reserves, New Phytol., 224, 625–631,, 2019. a, b

Dittmar, C., Fricke, W., and Elling, W.: Impact of Late Frost Events on Radial Growth of Common Beech (Fagus Sylvatica L.) in Southern Germany, Europ. J. Forest Res., 125, 249–259,, 2006. a, b, c, d

Nychka, D., Furrer, R., Paige, J., and Sain, S.: Fields: Tools for Spatial Data, (last access: 7 March 2023), 2021. a

Dulamsuren, C., Hauck, M., Kopp, G., Ruff, M., and Leuschner, C.: European Beech Responds to Climate Change with Growth Decline at Lower, and Growth Increase at Higher Elevations in the Center of Its Distribution Range (SW Germany), Trees, 31, 673–686,, 2017. a

Duveneck, M. J. and Thompson, J. R.: Climate Change Imposes Phenological Trade-Offs on Forest Net Primary Productivity, J. Geophys. Res.-Biogeo., 122, 2298–2313,, 2017. a

Ellenberg, H., Leuschner, C., and Dierschke, H.: Vegetation Mitteleuropas mit den Alpen: in ökologischer, dynamischer und historischer Sicht; 203 Tabellen, no. 8104 in UTB Botanik, Ökologie, Agrar- und Forstwissenschaften, Geographie, Verlag Eugen Ulmer, Stuttgart, 6., vollständig neu bearbeitete und stark erweiterte auflage, ISBN 978-3-8252-8104-5, ISBN 978-3-8001-2824-2, 2010. a

European Space Agency: Copernicus DEM,, 2024. a

Firke, S.: janitor: Simple tools for examining and cleaning dirty data, (last access: 3 March 2023), 2021. a

Gampe, D., Zscheischler, J., Reichstein, M., O'Sullivan, M., Smith, W. K., Sitch, S., and Buermann, W.: Increasing Impact of Warm Droughts on Northern Ecosystem Productivity over Recent Decades, Nat. Clim. Change, 11, 772–779,, 2021. a

Gazol, A., Camarero, J. J., Colangelo, M., de Luis, M., Martínez del Castillo, E., and Serra-Maluquer, X.: Summer Drought and Spring Frost, but Not Their Interaction, Constrain European Beech and Silver Fir Growth in Their Southern Distribution Limits, Agr. Forest Meteorol., 278, 107695,, 2019. a

Graves, S., Piepho, H.-P., and Dorai-Raj, L. S. with help from S.: multcompView: Visualizations of paired comparisons, (last access: 7 February 2020), 2019. a

Grossman, J. J.: Phenological Physiology: Seasonal Patterns of Plant Stress Tolerance in a Changing Climate, New Phytol., 237, 1508–1524,, 2023. a

Hickler, T., Smith, B., Sykes, M. T., Davis, M. B., Sugita, S., and Walker, K.: Using a Generalized Vegetation Model to Simulate Vegetation Dynamics in Northeastern USA, Ecology, 85, 519–530,, 2004. a

Hickler, T., Vohland, K., Feehan, J., Miller, P. A., Smith, B., Costa, L., Giesecke, T., Fronzek, S., Carter, T. R., Cramer, W., Kühn, I., and Sykes, M. T.: Projecting the Future Distribution of European Potential Natural Vegetation Zones with a Generalized, Tree Species-Based Dynamic Vegetation Model: Future Changes in European Vegetation Zones, Glob. Ecol. Biogeogr., 21, 50–63,, 2012. a, b, c

Hijmans, R. J.: terra: Spatial data analysis, (last access: 3 March 2023), 2022. a

Hufkens, K., Friedl, M. A., Keenan, T. F., Sonnentag, O., Bailey, A., O'Keefe, J., and Richardson, A. D.: Ecological Impacts of a Widespread Frost Event Following Early Spring Leaf-Out, Glob. Change Biol., 18, 2365–2377,, 2012. a, b, c, d, e

Inouye, D.: The Ecological and Evolutionary Significance of Frost in the Context of Climate Change, Ecol. Lett., 3, 457–463,, 2000. a, b

Jevšenak, J.: Daily Climate Data Reveal Stronger Climate-Growth Relationships for an Extended European Tree-Ring Network, Quaternary Sci. Rev., 221, 105868,, 2019. a

Keenan, T. F., Gray, J., Friedl, M. A., Toomey, M., Bohrer, G., Hollinger, D. Y., Munger, J. W., O'Keefe, J., Schmid, H. P., Wing, I. S., Yang, B., and Richardson, A. D.: Net Carbon Uptake Has Increased through Warming-Induced Changes in Temperate Forest Phenology, Nat. Clim. Change, 4, 598–604,, 2014. a

Kenk, G. and Guehne, S.: Management of Transformation in Central Europe, Forest Ecol. Manag., 151, 107–119,, 2001. a, b, c

Kollas, C., Körner, C., and Randin, C. F.: Spring Frost and Growing Season Length Co-Control the Cold Range Limits of Broad-Leaved Trees, J. Biogeogr., 41, 773–783,, 2014a. a, b

Kollas, C., Randin, C. F., Vitasse, Y., and Körner, C.: How Accurately Can Minimum Temperatures at the Cold Limits of Tree Species Be Extrapolated from Weather Station Data?, Agr. Forest Meteorol., 184, 257–266,, 2014b. a

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, J. Ecol., 104, 1076–1088,, 2016. a

Kramer, K., Ducousso, A., Gömöry, D., Hansen, J. K., Ionita, L., Liesebach, M., Lorenţ, A., Schüler, S., Sulkowska, M., de Vries, S., and von Wühlisch, G.: Chilling and Forcing Requirements for Foliage Bud Burst of European Beech ( Fagus Sylvatica L.) Differ between Provenances and Are Phenotypically Plastic, Agr. Forest Meteorol., 234/235, 172–181,, 2017. a, b, c

Lenz, A., Hoch, G., Vitasse, Y., and Körner, C.: European Deciduous Trees Exhibit Similar Safety Margins against Damage by Spring Freeze Events along Elevational Gradients, New Phytol., 200, 1166–1175,, 2013. a

Lenz, A., Hoch, G., Körner, C., and Vitasse, Y.: Convergence of Leaf-out towards Minimum Risk of Freezing Damage in Temperate Trees, Funct. Ecol., 30, 1480–1490,, 2016. a

Lindeskog, M., Arneth, A., Miller, P., Nord, J., Mischurov, M., Olin, S., Schurgers, G., Smith, B., Wårlind, D., and past LPJ-GUESS contributors: LPJ-GUESS Release v4.0.1 model code (4.0.1), Zenodo [code],, 2017. a, b

Lindeskog, M., Smith, B., Lagergren, F., Sycheva, E., Ficko, A., Pretzsch, H., and Rammig, A.: Accounting for forest management in the estimation of forest carbon balance using the dynamic vegetation model LPJ-GUESS (v4.0, r9710): implementation and evaluation of simulations for Europe, Geosci. Model Dev., 14, 6071–6112,, 2021. a

Liu, Q., Piao, S., Janssens, I. A., Fu, Y., Peng, S., Lian, X., Ciais, P., Myneni, R. B., Peñuelas, J., and Wang, T.: Extension of the Growing Season Increases Vegetation Exposure to Frost, Nat. Commun., 9, 426,, 2018. a

Lloret, F., Keeling, E. G., and Sala, A.: Components of Tree Resilience: Effects of Successive Low-Growth Episodes in Old Ponderosa Pine Forests, Oikos, 120, 1909–1920,, 2011. a, b

Ma, Q., Huang, J.-G., Hänninen, H., and Berninger, F.: Divergent Trends in the Risk of Spring Frost Damage to Trees in Europe with Recent Warming, Glob. Change Biol., 25, 351–360,, 2019. a, b, c, d

Massicotte, P. and South, A.: rnaturalearth: World map data from natural earth, Tech. Rep., (last access: 25 January 2024), 2023. a

Matsui, T., Eguchi, H., and Mori, K.: Control of Dew and Frost Formations on Leaf by Radiative Cooling, Environ. Con. Biol., 19, 51–57,, 1981. a

Medlyn, B. E., Duursma, R. A., and Zeppel, M. J. B.: Forest Productivity under Climate Change: A Checklist for Evaluating Model Studies, WIREs Climate Change, 2, 332–355,, 2011. a

Menzel, A., Estrella, N., and Fabian, P.: Spatial and Temporal Variability of the Phenological Seasons in Germany from 1951 to 1996, Glob. Change Biol., 7, 657–666,, 2001. a

Menzel, A., Helm, R., and Zang, C.: Patterns of Late Spring Frost Leaf Damage and Recovery in a European Beech (Fagus Sylvatica L.) Stand in South-Eastern Germany Based on Repeated Digital Photographs, Front. Plant Sci., 6, 1–13,, 2015. a, b, c, d

Meyer, B. F.: LPJ-GUESS Model code for “Frost matters: Incorporating late-spring frost in a dynamic vegetation model regulates regional productivity dynamics in European beech forests”, Zenodo [code],, 2024a. a

Meyer, B. F.: Reproducible analysis pipeline for “Frost matters: Incorporating late-spring frost in a dynamic vegetation model regulates regional productivity dynamics in European beech forests”, Zenodo [code],, 2024b. a

Meyer, B. F.: Data needed to reproduce analysis from “Frost matters: Incorporating late-spring frost in a dynamic vegetation model regulates regional productivity dynamics in European beech forests”, Zenodo [data set],, 2024c. a

Meyer, B. F., Buras, A., Rammig, A., and Zang, C. S.: Higher Susceptibility of Beech to Drought in Comparison to Oak, Dendrochronologia, 64, 125780,, 2020. a

Morin, X., Lechowicz, M. J., Augspurger, C., O'keefe, J., Viner, D., and Chuine, I.: Leaf Phenology in 22 North American Tree Species during the 21st Century, Glob. Change Biol., 15, 961–975,, 2009. a

Müller, K.: here: A simpler way to find your files, (last access: 3 March 2023), 2020. a

Neuner, G.: Frost Resistance in Alpine Woody Plants, Front. Plant Sci., 5, 654,, 2014. a

Nolè, A., Rita, A., Ferrara, A. M. S., and Borghetti, M.: Effects of a Large-Scale Late Spring Frost on a Beech (Fagus Sylvatica L.) Dominated Mediterranean Mountain Forest Derived from the Spatio-Temporal Variations of NDVI, Ann. Forest Sci., 75, 1–11,, 2018. a, b

Pebesma, E.: Simple Features for R: Standardized Support for Spatial Vector Data, The R Journal, 10, 439–446,, 2018. a

Pedersen, T. L.: patchwork: The composer of plots, (last access: 3 March 2023), 2020. a

Pedersen, T. L. and Crameri, F.: scico: Colour palettes based on the scientific colour-maps, (last access: 3 March 2023), 2022. a

Pierce, D.: ncdf4: Interface to unidata netCDF (version 4 or earlier) format data files, (last access: 3 March 2023), 2022. a

Pretzsch, H., Schütze, G., and Uhl, E.: Resistance of European Tree Species to Drought Stress in Mixed versus Pure Forests: Evidence of Stress Release by Inter-Specific Facilitation, Plant Biol., 15, 483–495,, 2013. a

Príncipe, A., van der Maaten, E., van der Maaten-Theunissen, M., Struwe, T., Wilmking, M., and Kreyling, J.: Low Resistance but High Resilience in Growth of a Major Deciduous Forest Tree (Fagus Sylvatica L.) in Response to Late Spring Frost in Southern Germany, Trees, 31, 743–751,, 2017. a, b, c, d

R Core Team: R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, (last access: 31 October 2023), 2022. a

Rammig, A., Jönsson, A., Hickler, T., Smith, B., Bärring, L., and Sykes, M.: Impacts of Changing Frost Regimes on Swedish Forests: Incorporating Cold Hardiness in a Regional Ecosystem Model, Ecol. Model., 221, 303–313,, 2010. a

Rubio-Cuadrado, Á., Gómez, C., Rodríguez-Calcerrada, J., Perea, R., Gordaliza, G. G., Camarero, J. J., Montes, F., and Gil, L.: Differential Response of Oak and Beech to Late Frost Damage: An Integrated Analysis from Organ to Forest, Agr. Forest Meteorol., 297, 108243,, 2021. a, b, c

Sangüesa-Barreda, G., Di Filippo, A., Piovesan, G., Rozas, V., Di Fiore, L., García-Hidalgo, M., García-Cervigón, A. I., Muñoz-Garachana, D., Baliva, M., and Olano, J. M.: Warmer Springs Have Increased the Frequency and Extension of Late-Frost Defoliations in Southern European Beech Forests, Sci. Total Environ., 775, 145860,, 2021. a

Scharnweber, T., Manthey, M., Criegee, C., Bauwe, A., Schröder, C., and Wilmking, M.: Drought Matters – Declining Precipitation Influences Growth of Fagus Sylvatica L. and Quercus Robur L. in North-Eastern Germany, Forest Ecol. Manag., 262, 947–961,, 2011. a

Schütz, J. P.: Close-to-Nature Silviculture: Is This Concept Compatible with Species Diversity?, Forestry, 72, 359–366,, 1999. a, b

Schwartz, M. D.: Assessing the Onset of Spring: A Climatological Perspective, Phys. Geogr., 14, 536–550,, 1993. a, b

Schwartz, M. D., Ahas, R., and Aasa, A.: Onset of Spring Starting Earlier across the Northern Hemisphere, Glob. Change Biol., 12, 343–351,, 2006. a

Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J. O., Levis, S., Lucht, W., Sykes, M. T., Thonicke, K., and Venevsky, S.: Evaluation of Ecosystem Dynamics, Plant Geography and Terrestrial Carbon Cycling in the LPJ Dynamic Global Vegetation Model, Glob. Change Biol., 9, 161–185,, 2003. a

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: Vegetation Dynamics in Ecosystem Models, Glob. Ecol. Biogeogr., 10, 621–637,, 2001. a, b

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,, 2014. a, b, c

South, A., Michael, S., and Massicotte, P.: rnaturalearthhires: High resolution world vector map data from natural earth used in rnaturalearth, Tech. Rep., 2024. a

Sullivan, P. F., Pattison, R. R., Brownlee, A. H., Cahoon, S. M. P., and Hollingsworth, T. N.: Effect of Tree-Ring Detrending Method on Apparent Growth Trends of Black and White Spruce in Interior Alaska, Environ. Res. Lett., 11, 114007,, 2016. a

Sykes, M. and Prentice, I.: Climate Change, Tree Species Distributions and Forest Dynamics: A Case Study in the Mixed Conifer/Northern Hardwoods Zone of Northern Europe, Climatic Change, 34, 161–177,, 1996. a

Urbanski, S., Barford, C., Wofsy, S., Kucharik, C., Pyle, E., Budney, J., McKain, K., Fitzjarrald, D., Czikowsky, M., and Munger, J. W.: Factors Controlling CO2 Exchange on Timescales from Hourly to Decadal at Harvard Forest, J. Geophys. Res., 112, G02020,, 2007. a, b

Vitasse, Y., Lenz, A., Hoch, G., and Körner, C.: Earlier leaf-out rather than difference in freezing resistance puts juvenile trees at greater risk of damage than adult trees, J. Ecol., 102, 981–988,, 2014. a

Vitra, A., Lenz, A., and Vitasse, Y.: Frost Hardening and Dehardening Potential in Temperate Trees from Winter to Budburst, New Phytol., 216, 113–123,, 2017. a

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., Takahashi, K., Vaughan, D., Wilke, C., Woo, K., and Yutani, H.: Welcome to the tidyverse, J. Open Sour. Softw., 4, 1686,, 2019. a

Wilmking, M., van der Maaten-Theunissen, M., van der Maaten, E., Scharnweber, T., Buras, A., Biermann, C., Gurskaya, M., Hallinger, M., Lange, J., Shetti, R., Smiljanic, M., and Trouillier, M.: Global Assessment of Relationships between Climate and Tree Growth, Glob. Change Biol., 26, 3212–3220,, 2020. a

Xu, K., Wang, X., Liang, P., An, H., Sun, H., Han, W., and Li, Q.: Tree-Ring Widths Are Good Proxies of Annual Variation in Forest Productivity in Temperate Forests, Sci. Rep., 7, 1945,, 2017. a, b

Yao, Y., Joetzjer, E., Ciais, P., Viovy, N., Cresto Aleina, F., Chave, J., Sack, L., Bartlett, M., Meir, P., Fisher, R., and Luyssaert, S.: Forest fluxes and mortality response to drought: model description (ORCHIDEE-CAN-NHA r7236) and evaluation at the Caxiuanã drought experiment, Geosci. Model Dev., 15, 7809–7833,, 2022. a

Yousefpour, R., Augustynczik, A. L. D., Reyer, C. P. O., Lasch-Born, P., Suckow, F., and Hanewinkel, M.: Realizing Mitigation Efficiency of European Commercial Forests by Climate Smart Forestry, Sci. Rep., 8, 345,, 2018. a

Zang, C., Hartl-Meier, C., Dittmar, C., Rothe, A., and Menzel, A.: Patterns of Drought Tolerance in Major European Temperate Forest Trees: Climatic Drivers and Levels of Variability, Glob. Change Biol., 20, 3767–3779,, 2014. a

Zeileis, A. and Grothendieck, G.: Zoo: S3 Infrastructure for Regular and Irregular Time Series, J. Stat. Softw., 14, 1–27,, 2005. a

Zimmermann, J., Hauck, M., Dulamsuren, C., and Leuschner, C.: Climate Warming-Related Growth Decline Affects Fagus Sylvatica, But Not Other Broad-Leaved Tree Species in Central European Mixed Forests, Ecosystems, 18, 560–572,, 2015. a

Zohner, C. M., Mo, L., Renner, S. S., Svenning, J.-C., Vitasse, Y., Benito, B. M., Ordonez, A., Baumgarten, F., Bastin, J.-F., Sebald, V., Reich, P. B., Liang, J., Nabuurs, G.-J., De-Miguel, S., Alberti, G., Antón-Fernández, C., Balazy, R., Brändli, U.-B., Chen, H. Y. H., Chisholm, C., Cienciala, E., Dayanandan, S., Fayle, T. M., Frizzera, L., Gianelle, D., Jagodzinski, A. M., Jaroszewicz, B., Jucker, T., Kepfer-Rojas, S., Khan, M. L., Kim, H. S., Korjus, H., Johannsen, V. K., Laarmann, D., Lang, M., Zawila-Niedzwiecki, T., Niklaus, P. A., Paquette, A., Pretzsch, H., Saikia, P., Schall, P., Šebeň, V., Svoboda, M., Tikhonova, E., Viana, H., Zhang, C., Zhao, X., and Crowther, T. W.: Late-Spring Frost Risk between 1959 and 2017 Decreased in North America but Increased in Europe and Asia, P. Natl. Acad. Sci. USA, 117, 12192–12200,, 2020. a

Short summary
Late-spring frost (LSF), critically low temperatures when trees have already flushed their leaves, results in freezing damage leaving trees with reduced ability to perform photosynthesis. Forests with a high proportion of susceptible species like European beech are particularly vulnerable. However, this process is rarely included in dynamic vegetation models (DVMs). We show that the effect on simulated productivity and biomass is substantial, warranting more widespread inclusion of LSF in DVMs.
Final-revised paper