Integration of a Frost Mortality Scheme Into the Demographic Vegetation Model FATES

Frost is damaging to plants when air temperature drops below their tolerance threshold. The set of mechanisms used by cold‐tolerant plants to withstand freezing is called “hardening” and typically take place in autumn to protect against winter damage. The recent incorporation of a hardening scheme in the demographic vegetation model FATES opens up the possibility to investigate frost mortality to vegetation. Previously, the hardening scheme was used to improve hydraulic processes in cold‐tolerant plants. In this study, we expand upon the existing hardening scheme by implementing hardiness‐dependent frost mortality into CLM5.0‐FATES to study the impacts of frost on vegetation in temperate and boreal sites from 1950 to 2015. Our results show that the original freezing mortality approach of FATES, where each plant type had a fixed freezing tolerance threshold—an approach common to many other dynamic vegetation models, was restricted to predicting plant type distribution. The main results emerging from the new scheme are a high autumn and spring frost mortality, especially at colder sites, and increasing mid‐winter frost mortality due to global warming, especially at warmer sites. We demonstrate that the new frost scheme is a major step forward in dynamically representing vegetation in ESMs by for the first time including a level of frost tolerance that is responding to the environment and includes some level of cost (implicitly) and benefit. By linking hardening and frost mortality in a land surface model, we open new ways to explore the impact of frost events in the context of global warming.

physiological and structural changes, such as the synthesis of antifreeze proteins and changes in lipid composition to increase membrane fluidity (Bansal et al., 2016;Chang et al., 2021). This enables plants to (a) tolerate extracellular formation of ice and the resulting cellular dehydration (Janská et al., 2010;Levitt, 1980), and, (b) avoid the formation of interstitial ice crystals by regulating nucleation (Wisniewski et al., 2014). The avoidance of interstitial ice crystallization can be achieved by keeping tissues isolated from the cold or by decreasing the temperature of ice nucleation. When anti-freeze proteins such as dehydrins are synthesized, the nucleation point of ice in tissues can be suppressed down to −38°C (Hanin et al., 2011). Then, at temperatures between −20°C and −30°C, the formation of intracellular glass, also named vitrification, further enables cold-acclimated woody plants to develop a resistance to much lower temperatures. Cold acclimation usually takes place during autumn, while its reversal occurs when temperatures start to rise again in spring to enable the reactivation of plant metabolism.
The maximum level to which plants can acclimate, as well as the timing and the rate at which they harden, is under genetic control and varies between and within species, depending on their developmental stage (Chang et al., 2021;Johansson et al., 2015;Sakai & Larcher, 1987). If moved to grow in the same outdoor location, like a common garden, northerly and highland ecotypes are likely to start hardening earlier than those originating from southerly or low latitude locations (Stevenson, 1994). After exposure to moderately warm conditions, a decrease in frost hardiness takes place (Pagter & Arora, 2013;Vyse et al., 2019). The rate and temperature range at which plants deharden or deacclimate depends on the development status and genotype of species (Arora & Rowland, 2011).
The process of deacclimation is typically faster than acclimation (Arora et al., 1992). An extreme example shows that leaves of Solanum commersonii (wild potato species) only required 1 day to completely lose their freezing tolerance when exposed to 20°C (Chen & Li, 1980). In comparison, if exposed to 2°C, 15 days were needed. Observations of Pinus sylvestris indicate that hardening rates are around −1°C per day, while dehardening rates are ca. 2°C per day (Beck et al., 2004;Repo, 1991).
Frost mortality typically occurs at northern latitudes or high altitudes, following warm pulses and abrupt transitions from spring or even summer-like temperatures and back to typical winter conditions. During such events, cold tolerance can be dramatically disturbed. For example, in northern Scandinavia, temperatures rose up to 7°C in December 2007, resulting in loss of snow cover and exposure of vegetation to first warm and then returning cold temperatures (Bokhorst et al., 2009). Such "warm spell" events generate spring-like conditions which in extreme cases trigger a complete separation from a plant's freezing tolerance and further expose vegetation by melting of the protective snow cover (Bokhorst et al., 2011). Damage to trees and shrubs are likely to occur if temperature abruptly returns to freezing following a warm spell. Plants exposed to temperatures below their tolerance threshold will suffer from (a) mechanical stress, due to the intercellular ice crystallization and fragility of the tissues, and (b) osmotic dehydration of cells due to the freezing of intercellular water (Janská et al., 2010;Levitt, 1980). If vegetation is damaged during several successive winters, shifts in ecosystem composition may be observed (Zhao et al., 2017).
Global warming is changing the boundary conditions that determine the winter survival of plants (Saxe et al., 2001). As the Earth warms, winters in cold regions become milder, growing seasons longer, and the snowpack thinner. However, as temperatures rise and become more variable, the frequency and severity of extreme weather events are expected to increase, including the occurrence of winter warm spells (Vikhamar-Schuler et al., 2016). Research involving aspects of freeze tolerance and winter survival of plants increases in relevance in times of rapid climate change, such that improved understanding of frost mortality is of practical importance today.
Global land surface models are the terrestrial elements of Earth System Models, which are our primary tools for predicting climate change, the responses of the biosphere to change, and the feedbacks from land surface processes onto the climate system. They provide mechanistic linkages between systems that interact to control these processes-hydrology, micrometeorology, biogeochemistry, vegetation physiology, ecology and cryospheric processes . Frost mortality is often poorly represented or totally absent in land surface models. In the more advanced dynamic vegetation models, climate envelopes are typically represented by assigning a minimum temperature (between 2.5°C and −80°C) to which a plant functional type (PFT) is tolerant and can survive. This is also the case for the FATES module of the Community Land Model (CLM). FATES is a size-and age-structured representation of vegetation dynamics which can be coupled to a land surface model or an Earth System Model-ESM (Fisher et al., 2015;Koven et al., 2020). The temperature thresholds used in the current freezing mortality parameterization of FATES were measured by Sakai and Weiser (1973) and Albani et al. (2006). Sakai and Weiser (1973) concluded that habitat is the main driver of plants' cold tolerance. They showed that conifer and deciduous trees of the boreal region can resist temperatures as low as −80°C, while in the Pacific and south-eastern coastal regions of North America, vegetation only tolerates temperatures down to ca. −15°C. In the default FATES model, it uses these static temperature thresholds to predict freezing mortality. Plants with a limited freezing tolerance (i.e., 2.5°C) are unlikely to grow in temperate and boreal regions, while plants with a strong freezing tolerance (i.e., −80°C) may never incur injury. In reality, the ability to withstand frost varies during the course of a year as plants slowly cold-acclimate during the onset of winter and de-acclimate in spring or when winter warm spells interrupt dormancy (Leinonen et al., 1997;Rammig et al., 2010).
The recent incorporation of a hardening scheme into the CLM5-FATES model, to better represent winter plant hydraulics (Lambert et al., 2022) (hydro-hardening scheme) opens up the possibility to include the hardiness level of a PFT into the calculation of frost mortality. In the "hydro-" hardening scheme parameterization (Lambert et al., 2022), the hardiness level is (a) variable throughout the course of a year, (b) is similar for all PFTs and (c) depends on the climate of its location and the seasonal changes in solar radiation. The "derived" hardiness level is then used to "influence" conductances along the soil-plant-atmosphere continuum and the rate for hydraulic failure mortality ( Figure S10 in Supporting Information S1). In this study we expand on the hardening calculation of the hydro-hardening scheme, by modifying it to also become a frost-hardening scheme, from which we calculate freezing mortality ( Figure S10 in Supporting Information S1).
While the hardening scheme is common to both hydro-as well as frost-hardening, we have here turned off the hydraulics functionality in FATES to focus on the hardening scheme's implications for an updated freezing mortality. We further develop the hardening scheme by using the freezing tolerance threshold (of the default FATES) to predict the maximum hardiness level of a PFT. We then implement a new frost scheme by replacing the fixed freezing tolerance threshold in the original freezing mortality calculation of FATES (Albani et al., 2006) with the time-and now also PFT-dependent hardiness level.
In addition, we analyze the new frost scheme to assess its impact on vegetation dynamics. While we assume that the use of a static PFT-dependent temperature threshold might be a good approximation to predict distribution across climate gradients, we hypothesize that the current freezing mortality parameterization is not capable of accurately identifying the impact, and variability, of damaging frost events nor the complex temporal dynamics of the risk of frost mortality. We further hypothesize that the new frost scheme (a) is necessary to model (realistic) frost mortality in northern regions of the globe, and (b) may lead to significant changes of PFT distribution by enhancing the competitiveness of frost-tolerant PFTs, leading to more reliable projections of how northern ecosystems respond to climate change.  (Golaz et al., 2019) and Norwegian ESM (Seland et al., 2020). CLM has expanded and seen considerable developments over the last decades to provide a comprehensive platform for researchers to study weather, climate change, and climate variability in response to terrestrial processes (Lawrence et al., 2019). CLM represents a wide range of terrestrial biogeophysical and biogeochemical processes, including the hydrological cycle and snow features, ecosystem dynamics, river transport, energy and greenhouse gas fluxes, and land use change, among others.

Materials and Methods
FATES is a vegetation demographic model, which is coupled to the CLM to better represent vegetation processes and dynamics. It builds on the disturbance dynamics of the Ecosystem Demography model of Moorcroft et al. (2001), wherein individual trees are grouped into cohorts based on their size and PFT and divided into successional stages. The model includes allometry and allocation schemes, flexible PFT trait definition, logging, and plant hydrodynamics among others. A detailed description can be found in Fisher et al. (2010Fisher et al. ( , 2015, Holm et al. (2020), and Koven et al. (2020). FATES incorporates a large set of plant mortality mechanisms, including background, freezing, carbon starvation, hydraulic failure, disturbance impact, logging, and fire mortalities.
In the default version of FATES, a PFT experiences freezing mortality when the daily mean vegetation temperature (Tveg) within a grid-cell drops below the freezing tolerance temperature threshold (PFTfreezetol) of that same PFT (Table 1). The rate at which a PFT dies is proportional to the difference between the vegetation temperature and the freezing tolerance threshold. In a cohort-based model, mortality is a rate predicated on a predictor variable, not a discrete event-as it might be in an individual-based model or in the real world. Thus, a scaling coefficient is required to generate cohort-wide loss of individuals from the freezing status predictor (Fisher et al., 2010). For freezing mortality, the delta is multiplied by a fixed scalar (PFTcoldstress) which is set to 3 (fraction per year) for all PFTs, which effectively means ∼0.8% per day (these values can be considered as a tuning factor in the use of all cohort-based vegetation schemes). Finally, a 5°C buffer (FrostBuffer) is used, so that mortality starts at 0% of individuals per day when vegetation temperature is 5°C above the tolerance threshold, and linearly increases to ∼0.8% of individuals per day as vegetation temperature approaches the freezing tolerance threshold (Equation 1, Figure S9 in Supporting Information S1). Whenever vegetation temperature is equal or below the tolerance threshold, mortality remains at ∼0.8% of individuals per day. More informa tion about the parameterization of the density-independent freezing mortality can be found in Albani et al. (2006), and the measurements of threshold temperatures for cold tolerance are reported in Sakai and Weiser (1973).

Hardening Scheme
The hardening scheme implemented in FATES by Lambert et al. (2022), was based on the work by Rammig et al. (2010). The scheme by Rammig et al. (2010) was developed to fit the climate of central Sweden (Farstanäs) and measurements on Norway spruce (Picea abies). During its implementation into FATES, the hardening scheme was modified to function globally. However, it is still lacking the ability to differentiate PFTs. In the current hardening scheme, the hardiness level (identical for all PFTs) is updated on a daily basis, depending on the mean daily temperature (T). For each model day, the hardiness level (HD) is calculated using three functions: the target hardiness (TH), the hardening rate (HR), and the dehardening rate (DR) (Equations 3-5, Figure S8 in Supporting Information S1). The first function is the target hardiness and is used to define whether the hardiness level will increase or decrease during the current time-step. If the target hardiness is lower than the hardiness of the previous day (HDP), the hardening rate is subtracted from the hardiness level of the previous day, meaning that the hardiness level decreases. By contrast, if the target hardiness is higher than the hardiness level of the previous day, the dehardening rate is added to the hardiness level of the previous day, and hence, the hardiness level increases during the current time-step (Equation 2).
The scheme incorporates a site-specific temperature index T5, defined as the 5-year running mean of the yearly minimum daily mean temperature. This is necessary to make the scheme applicable at a global scale. T5 governs the maximum hardiness level (H MAX ) and thereby enables vegetation of colder sites to harden down to a lower threshold than warmer locations (Equation 6). H MAX is also used in the calculation of HR and DR, so that plants of colder locations have the capacity to harden and deharden faster than plants growing in warmer locations. This is necessary to enable plants of cold locations to quickly reach a sufficiently low hardiness level and to prevent them to be systematically killed in autumn.  (Sakai & Weiser, 1973) 1. The target hardiness: 2. The hardening rate 3. The dehardening rate In Rammig et al. (2010), hardening was prevented until the 210th Julian day and a fixed growing degree day threshold was reached. In the implementation by Lambert et al. (2022), another location-specific adaptation was to introduce a threshold (DaylThresh, in seconds) to determine the onset of the hardening period (Equation 7 and Figure 1). After the daylength threshold is crossed for the first time after the summer solstice, dehardening is not possible anymore until days start lengthening again. This addition enables high latitude sites to start hardening when the daylength is longer than at lower latitudes. To avoid abnormally large mortality rates during autumn at colder sites of a given latitude, DaylThresh was made dependent on the temperature index T5 of the corresponding site. It is well-known that the decrease in temperature (mainly) and the decrease in daylength (to a lesser extent) (Aronsson, 1975;Beck et al., 2004) are crucial to determine the onset of dormancy and hardening in overwintering plants. DaylThresh takes both these environmental cues into account.

New Frost-Hardening Mortality Scheme
The previous implementation of a hardening scheme into FATES by Lambert et al. (2022) aimed to improve the representation of hydraulic processes in the Arctic and boreal region during winter (hydro-hardening). However, the hardening scheme by Rammig et al. (2010) was initially meant to be coupled to a frost damage model. The frost damage model by Rammig et al. (2010) uses the hardiness level as an input and yields a growth-reducing factor whenever the daily minimum temperature goes below the hardiness level of that day. A large spread between the minimum temperature and HD indicates a stronger damaging frost event which translates into a higher growth reducing factor.
In this study, we take advantage of the hardening scheme in FATES to update the current parameterization of freezing mortality by Albani et al. (2006). To achieve this, we first expand the existing hardening scheme (Lambert et al., 2022) by considering the tolerance threshold of PFTs (PFTFreezetol) in the calculation of the hardiness level (Equation 8). By doing so, each PFT is associated with a unique hardiness level depending on its tolerance to cold. We then combined the updated version of the hardening scheme with the approach by Albani et al. (2006) to calculate a new frost mortality, in the sense that we incorporate the time-dependent (as in Rammig et al. (2010)) and now also PFT-dependent hardiness level into the mortality calculation. We call the updated version of the hardening scheme combined with the new approach to calculate freezing mortality the "frost-hardening scheme" ( Figure S10 in Supporting Information S1). In the original FATES, PFTfreezetol was directly used to predict freezing mortality, but now it only represents the maximum possible value for H MAX of the corresponding PFT (Equation 8). This means that in the new frost-hardening scheme, frost mortality becomes a function of HD.
The new implementation of frost-hardening constrains the hardiness level between −2°C and −70°C as the most extreme values and allows it to vary within that range. The default H MAX is defined as T5-10°C (Equation 8), which means that the maximum hardiness level is 10°C below the 5-year running mean of the yearly minimum daily mean temperature. We also tested two other H MAX thresholds: T5-5°C and T5-15°C. MAX = min(max(PFTFreezetol, max( 5, −60) − 10), MIN) Equation 9 forms the core of the frost scheme, where FRmort is the actual frost mortality in fraction of individuals per day. As in the original freezing mortality approach of FATES (Equation 1), FrostBuffer is set to 5°C and PFTcoldstress scalar equals ∼0.008 (fraction per day). As in Rammig et al. (2010), the minimum daily atmospheric temperature (T MIN ), and the hardiness level (HD) are used to predict mortality.

Model Experiments
Model experiments were carried out using CLM5.0-FATES, with fully prognostic variables for vegetation, litter, and soil. The atmospheric forcing to drive model simulations was extracted from the Global Soil Wetness Project 3 version 1 (GSWP3v1) data set (Dirmeyer et al., 2006). GSWP3v1 is a 3 hourly resolution observed climate data set on a global 0.5° × 0.5° grid. The data set covers the period 1901 to 2015, which makes it a valuable data set to look into the evolution of frost events throughout the 21st Century and beyond. For our study, we used near-surface specific humidity, daily mean temperature, near-surface wind magnitude, downwelling longwave radiation, total precipitation, shortwave downwelling radiation, and surface air pressure from 1921 to 2015.
We show the functioning of the frost-hardening mortality scheme for single site simulations at Soběšice (in the Czech Republic), Farstanäs (in Sweden), and Spasskaya Pad (in Russia). These sites were selected to illustrate the performance of the frost-hardening scheme in different climates (Figure 2). For each site, we tested the behavior of the frost-hardening scheme for different PFTs with various PFTfreezetol parameters. We decided to run single PFT simulations as calibration of cohort-based vegetation demographic models to represent coexistence remains challenging (Buotte et al., 2021;Koven, 2019), and this is an orthogonal problem to the present study. We decreased the value of PFTfreezetol for broadleaf cold-deciduous extratropical trees from −30°C to −80°C in this study as it is the only PFT representing broadleaf deciduous trees in the boreal regions, and survives much lower temperatures than −30°C. We selected the following PFTs to evaluate the new frost mortality scheme across a broad range of freezing tolerances: The occurrence and magnitude of damaging frost events was studied from 1950 to 2015. The frost mortality variable described in this study represents the percentage of individuals dying per day, and is therefore independent of stem density and successional competitive dynamics. In other words, the "spinup" period , during which vegetation tends toward equilibrium, does not influence the frost mortality rates shown in this paper.

General Functioning of the Frost Scheme
In the control version of FATES, a PFT with a freezing tolerance threshold of −80°C (such as broadleaf deciduous extratropical trees) does not experience a single damaging frost event between 2010 and 2013, at any of the three high latitude sites shown in Figure 3. However, while using the hardiness scheme to predict frost mortality for the same PFT, several events take place during this 3-year period.
At Soběšice, it is not uncommon for winter temperatures to occasionally increase beyond 2.5°C, which is the minimum required value for dehardening. It follows that the hardiness level oscillates throughout the cold season, which makes it possible that it is crossed in the middle of the winter, leading to repeated mid-winter damaging frost events. At Spasskaya Pad, in contrast, temperatures generally remain well below 0°C during winter, which means that the hardiness level does not increase until spring. While the occurrence of damaging frost events is almost non-existent in mid-winter, they are possible during the shoulder seasons when plants are not fully hardened. Farstanäs is characterized by colder winters than Soběšice, which means that mid-winter dehardening and mid-winter mortality from frost events is much rarer. On the other hand, Farstanäs is warmer than Spasskaya Pad and the cold seasons are usually much shorter. This results in reduced spring and autumn frost events in Farstanäs.

Applicability of the Scheme Across Biomes
The functioning of the frost scheme across the globe was tested by running CLM5.0-FATES at multiple locations in different parts of the world (Figure 2). Figure 4 illustrates how the hardiness model behaves in the three main locations of this study (Figures 4a-4c), as well as sanity checks in a tropical location where the daily temperature rarely drops to levels that trigger hardening (Caxiuanã, Figure 4d), at a location in the southern hemisphere with an oceanic climate (Temuco, Figure 4e), and lastly, a site at a lower latitude than the three main sites with a continental climate and cold winters (Harvard Forest, Figure 4f). Several PFTs with different freezing tolerance thresholds were modeled for each study site. Although hardiness levels of the study areas partly overlap, we show that the hardiness level, and hence the frost mortality parameterization, keeps its function as a distribution predictor in the new frost scheme. For example, BET trees do not have the capacity to acclimate to cold weather and will incur mortality as soon as minimum temperature falls below −2°C. This makes it challenging for these trees to compete for water and nutrients at locations such as Soběšice and impossible to survive in Spasskaya Pad. However, in contrast to the current freezing mortality parameterization, PFTs can experience frost during the shoulder and summer seasons since the hardiness level, even of the most tolerant PFTs, starts at the minimum hardiness level of −2°C.

Sensitivity to the Maximum Hardiness Level
Frost mortality rates are sensitive to the definition of the maximum hardiness level (the degree to which plants are capable of hardening). H MAX is unique to each PFT as it depends on PFTfreezetol and varies from site to site as a function of the 5-year running mean of temperature.
Given that this feature of the model is at present unconstrained by physiological observations, we here tested two additional methods for defining H MAX to illustrate the sensitivity of the model to this parameter, where either 5°C or 15°C was subtracted from T5 to define H MAX ( Figure 5 and Figure S1 in Supporting Information S1). It is clear that a higher H MAX (T5-5°C) systematically results in higher mortality rates. The difference between T5-5°C and T5-10°C is much larger than between T5-10°C and T5-15°C. Additionally, at the two warmest sites (Figures 5a and 5b), the spread between T5-5°C and the other thresholds is considerably larger than in Spasskaya Pad. Figure 5 also shows that a PFT with a freezing tolerance of −80°C will experience more frost at Spasskaya Pad than at Soběšice, and more at Soběšice than at Farstanäs, similar to Figure 3.

Monthly Distribution of Frost Events
To illustrate the annual distribution of frost events, Figure 6 shows the 30-year average  of simulated frost mortality rates per month for BDE (broadleaf deciduous extratropical) trees. BEE (broadleaf evergreen extratropical) and NEE (needleleaf evergreen extratropical) trees are shown in Supporting Information S1 ( Figures S2 and S3). Mortality rates were first summed over each month of the year, then averaged over the 30-year period. The total number of days with frost mortality (>0% of individual day −1 ) during the 30-year period is indicated above each bar. In the control FATES simulation, (a) BDE (broadleaf deciduous extratropical) trees (−80°C) do not experience mortality from a single frost event during 30 years at any of our sites, including Siberia ( Figure 6), (b) BEE (broadleaf evergreen extratropical) trees (−30°C) quickly die out at Spasskaya Pad, but never experience frost mortality at the other two sites ( Figure S2 in Supporting Information S1), and (c) only NEE (needleleaf evergreen extratropical) trees, for which the freezing tolerance (−55°C) is slightly above the minimum winter temperatures at Spasskaya Pad, the model simulates small amounts of frost mortality, but these are far from causing complete dieback of the PFT ( Figure S3c in Supporting Information S1).
In the new frost-hardening scheme, PFT distribution is still mainly predicted by PFTfreezetol since non-tolerant PFTs do not withstand the climate at Spasskaya Pad ( Figure S2c in Supporting Information S1). However, even the most tolerant PFTs can now experience frost mortality while they are in a dehardened state.
At Soběšice, most frost mortality occurs in the middle of the winter (December and January), and the amounts are the same for all PFTs since they have the same hardiness levels (Figure 6a, Figures S2a and S3a in Supporting Information S1). There are also several frost events during the shoulder seasons at Soběšice, especially in April and October. Of these three sites, Farstanäs is the site with the smallest number of damaging frost events, but December and January have some frost mortality (Figure 6b, Figures S2b, S3a, and S3b in Supporting Information S1).
In contrast to the two other sites, PFTs at Spasskaya Pad experience contrasting amounts of frost as two of them do not tolerate the low temperatures. The two least tolerant PFTs undergo near-complete winter mortality ( Figure S3c in Supporting Information S1), even causing BEE trees to die out ( Figure S6c in Supporting Information S1). This is why the monthly distribution of frost events is not shown for BEE trees ( Figure S2 in Supporting Information S1). For BDE trees (−80°C), frost mortality occurs most frequently in the shoulder seasons ( Figure 6c). In contrast, due to the strong hardiness levels in the middle of winter, no damaging frost event takes place from December to April (Figure 6c). These findings are corroborated by a study of Horvath et al. (2021) showing that spring and autumn time temperatures are critical for describing birch and pine distribution.

Trends of Frost Events
The trend in damaging frost events was assessed by comparing the first 30 years  of our model runs to the last 30 years   (Figure 7, Figures S4 and S5 in Supporting Information S1). In Figure S4 in Supporting Information S1, only results for Soběšice and Farstanäs are shown, since BEE (broadleaf evergreen extratropical) trees did not survive at Spasskaya Pad. Spasskaya Pad is the site with the highest frost mortality, which makes the results more robust (Figure 7c). The comparison of the early and late time periods shows that damaging frost events decrease strongly, especially during the warmest months of the year, and particularly in Spasskaya Pad where early summer frosts events are generally simulated to be frequent (Figure 7 and Figure S5 in Supporting Information S1). In Soběšice and Farstanäs, frost mortality in mid-winter months-especially January and December-increases noticeably (Figures 7a and 7b), underlining the complex dynamics of these processes in a warming climate.

PFT Living Biomass
In most cases (sites and PFTs), the application of the new frost-hardening scheme yielded similar PFT biomass to the biomass in the control FATES simulation (Figures S6a, S6b, S7a, and S7b in Supporting Information S1;  Figures 8a and 8b). The largest differences in vegetation biomass were seen at Spasskaya Pad ( Figure S7c in Supporting Information S1 and Figure 8c).

Conclusions
The original frost mortality implementation in FATES (Albani et al., 2006) uses hard-coded temperature thresholds to constrain the distribution of vegetation across biomes. These mortality thresholds prevent non-tolerant PFTs from becoming too competitive in northern climates, and in multi-PFT model simulations, to replace boreal PFTs. Despite the challenges to define those thresholds and the need to replace them by more process-based descriptions, this static approach is used in most dynamic vegetation models to define climate envelopes (Beigaitė  , 2022). Due to the presence of fixed temperature thresholds, this approach does not account for climate variations and the resilience of vegetation to frost is independent of weather conditions. In reality, variations in climate are major determinants of vegetation distribution across climate envelopes (Adams, 2009). As the world is currently facing fast climate changes, modifications of the large-scale dynamics in vegetation distribution are expected (Pearson et al., 2013).
In this study, we introduce an alternative scheme that fulfills the same role of governing PFT distribution, but at the same time introduces more dynamic physiologically based mechanisms. This new frost-hardening scheme has the capacity to simulate frost mortality when plants are unprepared to withstand freezing. The freezing tolerance threshold is still unique to each PFT but is also variable in time depending on the environmental cues of a given location. This is of crucial importance, as it enables the detection and prediction of the damaging effects of extreme winter events and late spring frost events.
We show that at some locations (Soběšice and Farstanäs), and for some PFTs (BDE, NEE and NDE trees), the new frost-hardening scheme will cause minor changes to total PFT biomass, while at other sites (e.g., Spasskaya Pad), much lower amounts of biomass were simulated during the spin-up (e.g., Figure S8c in Supporting Information S1). We argue that such situations may have large impacts on competition in multiple PFT simulations, but more work with parameter estimation of high-latitude PFTs in FATES remains necessary before it can be run realistically in competition mode (Buotte et al., 2021).
In the new frost-hardening model, the hardiness level of all PFTs ranges from −2 to −70°C. Mortality increases linearly from when minimum temperatures cross the hardiness level until they reach 5°C below the hardiness level. In the original freezing mortality approach, mortality would start 5°C above the PFTfreezetol parameter (i.e., between 2.5 and −80°C). While a PFT tolerant to 2.5°C would start experiencing freezing mortality at 7.5°C in the old scheme, the same PFT will in the new scheme start experiencing frost mortality at −2°C.
In the frost-hardening scheme, the temperature threshold at which a PFT experiences frost mortality is no longer fixed. Rather, frost mortality is triggered by the hardiness level which evolves dynamically as a function of the climatic conditions and involves a set of variables and parameters in its calculation, including the minimum and maximum hardiness level, the target hardiness function, hardening response rate, temperature range, rates of dehardening, as well as the daylength threshold before which no hardening is possible in autumn.
The sensitivity of the scheme to the hardening rate, dehardening rate, maximum hardiness level, and start of autumn (equivalent to the daylength threshold) were already discussed in the first version of the hardening scheme (Rammig et al., 2010). In short: they concluded that the hardening model was sensitive to variations in the hardening rate, and not sensitive to the dehardening rate and changes in the target hardiness. The frequency and magnitude of frost mortality decreased with an increasing hardening rate and vice-versa. While the overarching concept remains the same, our implementation is a considerable modification of the original hardening scheme (Rammig et al., 2010) to make it applicable globally and for multiple PFTs.
Unless a PFT is intolerant to the temperatures of its site, damaging frost events are short lived in the new scheme. They rarely last more than a few days, and their impact is limited to 0.8% of individuals per day. Therefore, the potential of impact of frost on vegetation is small in comparison to the parameterization of carbon starvation and hydraulic failure mortalities, which typically remain at high levels for several months during winter in northern regions (Lambert et al., 2022). We did not modify the amplitude parameter due to a lack of observations to evaluate against, but it can be increased when necessary. We have tried to compare our model to data sets on frost mortality based on tree ring data (Gurskaya, 2021), but these show damage to the xylem, not needle loss or mortality of saplings (which do not survive to be included in tree ring observations). This makes these data sets unsuitable for validation of our model. Observations of frost damage from forestry databases (like the Norwegian skogskader.no) are extremely sparse, and the data set is likely to contain multiple false negatives meaning frost damage events are underreported. We call for more quantitative and systematic observations on frost mortality (such as crown defoliation) to achieve a precise parameterization of the frost scheme. Rammig et al. (2010) validated their model by comparing average mortality rates across 122 sites in Sweden (mostly in southern and central Sweden, i.e., below 61°N). When comparing the performance of our adaptation of this scheme, we broadly find similar mortality rates across southern and central Sweden (1.5%, Table S1 in Supporting Information S1) as Rammig et al. (2.2%) and their observation estimate (1.1%) for the period 1999 to 2005. Over large areas, we therefore believe that our scheme provides plausible frost mortality rates. We show that our new frost scheme gives moderate amounts of frost mortality and is conservative at the studied sites. It is also an improvement upon the existing scheme in FATES which showed no mortality to PFTs in areas where frost damage is known to occur.
Research on cost and benefits of hardening is scarce and it has been suggested that there may not be any cost in the absence of freezing stress (Zhen et al., 2011). However, when freezing occurs, damage to membranes and photosynthetic enzymes and proteins may be costly to plants (Chang et al., 2021;Gray et al., 1997). When the hydraulic scheme of FATES and the hardening scheme are simultaneously activated, hardening becomes costly to plants as it reduces the conductivity between plant water pools and soil (Lambert et al., 2022). In FATES-Hydro, the soil-plant continuum is discretized in a series of water storage compartments with variable heights, volumes, water retention, and conducting properties (Christoffersen et al., 2016). Hence, in the Hydro scheme, water potentials are predicted for each plant organ (roots, stem, leaves), and fluxes between organs depend on conductance and water potential gradients. Lambert et al. (2022), extensively discuss how the decrease in conductivity between plant tissues mimics the effect of hardening on the decrease of xylem water fluxes (Gusta et al., 2005;Smit-Spinks et al., 1984;Steponkus, 1984). By diminishing root water exchanges (influx and efflux), photosynthesis rates remain low throughout the winter preventing high rates of productivity. Although these secondary effects might not have been intended in the first place, they seem to fit the behavior of hardy plants in the field, since they are most often also dormant, which means that their metabolism is slowed down (Chang et al., 2021;Havranek & Tranquillini, 1995). In contrast, in FATES without the Hydro scheme, plants are considered as a simple resistance to water and the soil water potential predicts plant water fluxes. In this study, the hydraulic module of FATES was turned off to focus on the frost mortality. As Hydro was turned off, so were the hydraulic constraints implemented by Lambert et al. (2022), meaning that hardening had no cost in the results presented in this paper.
The frost-hardening scheme frees new ways to explore the evolution of frost occurrences and intensities under changing climatic conditions. However, further research is necessary to better assess the distribution and impact of frost events on vegetation across the globe. Our results directly depend on the forcing, and therefore, we recommend further testing of several climate sources. Most importantly, we encourage further research to evaluate the hardening-frost scheme against field observations to better assess the magnitude and timing of frost events, which is why we designed the model framework to be flexible and easily adaptable. As winter weather is expected to become more variable (Cohen et al., 2018;Graham et al., 2017), increasing the frequency of mid-winter warm spells and frost events during the shoulder seasons, our hardening-based frost mortality scheme is a necessary model improvement to better quantify the impact of frost mortality on vegetation dynamics and potential species shifts with future climate change.
Data Availability Statement also like to thank the World Climate Research Programme for the GSWP3v1 reanalysis product. Finally, we would like to acknowledge our colleagues at NCAR and Berkeley lab for their precious help during model setup and manuscript preparation. RF acknowledges the support of the US Dept of Energy "Next Generation Ecosystem Experiment in the Tropics" project, as the EUH2020 4C project. The simulations were performed on the SAGA supercomputer operated by Sigma2 (project number NN2806K).