Exposure to cold temperature affects the spring phenology of Alaskan deciduous vegetation types

Temperature is a dominant factor driving arctic and boreal ecosystem phenology, including leaf budburst and gross primary production (GPP) onset in Alaskan spring. Previous studies hypothesized that both accumulated growing degree day (GDD) and cold temperature (chilling) exposure are important to leaf budburst. We test this hypothesis by combining both satellite and aircraft vegetation measurements with the Community Land Model Version 4.5 (CLM), in which the end of plant dormancy depends on thermal conditions (i.e. GDD). We study the sensitivity of GPP onset of different Alaskan deciduous vegetation types to a GDD model with chilling requirement (GC model) included. The default CLM simulations have a 1–12 d earlier day of year GPP onset over Alaska vegetated regions compared to satellite constrained estimates from the Polar Vegetation Photosynthesis and Respiration Model. Integrating a GC model into CLM shifts the phase and amplitude of GPP. During 2007–2016, mean GPP onset is postponed by 5 ± 7, 4 ± 8, and 1 ± 6 d over Alaskan northern tundra, shrub, and forest, respectively. The GC model has the greatest impact during warm springs, which is critical for predicting phenology response to future warming. Overall, spring GPP high bias is reduced by 10%. Thus, including chilling requirement in thermal forcing models improves northern high-latitude phenology, but leads to other impacts during the growing season which require further investigation.


Introduction
Vegetation phenology is crucial to the dynamic response of terrestrial ecosystems to climate change (Euskirchen et al 2014), since vegetation affects surface radiation, temperature, energy exchange, the hydrological cycle, and the carbon cycle (Myneni et al 1997, Schaefer et al 2005, Euskirchen et al 2006, Piao et al 2007, Bonan 2008, Sitch et al 2008, Jeong et al 2011, Euskirchen et al 2014. A better understanding and prediction of ecosystem phenology is essential for reducing uncertainties in modeling carbon, water, and energy cycles, and their feedback to climate (Levis and Bonan 2004). The model comparison of the North American Carbon Program model intercomparison suggests that compared to observations, most terrestrial biosphere models (TBMs) simulate earlier starts (i.e. typically two weeks) of the growing season in deciduous forests, and the predicted leaf onset dates under different climate scenarios are highly uncertain as a result of the dependence of predicted response to warming on model structures (Richardson et al 2012).
Further, early starts of the growing season could be the reason for the over-estimation of gross ecosystem photosynthesis during the spring transition in TBMs (Richardson et al 2012). Thus, accurate simulation of vegetation phenology remains a major challenge in TBMs.
In most TBMs, leaf onset of deciduous types could be triggered in different ways: temperature threshold (El Maayar et al 2002), growing degree day (GDD; Sitch et al 2003), GDD and radiation accumulation (Thornton et al 2002), and GDD and chilling requirement (Krinner et al 2005). Here, GDD refers to the accumulated mean surface air temperature above a given threshold (e.g. 0°C). Studies have shown that some temperate and boreal deciduous tree species need the exposure to certain ranges of cold temperature (i.e. chilling requirement) to break dormancy (Baldocchi and Wong 2008, Schwartz and Hanes 2010, Delpierre et al 2016. Here, the chilling requirement is species-specific (Laube et al 2014); for example, the number of chill days (with daily mean surface air temperature less than 5°C) before the break of dormancy is ∼100 for temperate deciduous trees (Jeong et al 2012). Modeling based studies also support that the chilling requirement is essential to deciduous species for the break of dormancy, regulating the timing of leaf budburst (White et al 1997, Krinner et al 2005, Jeong et al 2012, Richardson et al 2012, Clark et al 2014. Thus, there have been several efforts to move beyond a simple thermal forcing model approach (GDD models;White et al 1997, Sitch et al 2003 and test more sophisticated GDD models that account for chilling period (herein GC model; Krinner et al 2005, Caffarra et al 2011, Jeong et al 2012. For example, Jeong et al (2012) find that using a chilling model shifts the zero-crossing date for net carbon uptake by ∼11 d in temperate deciduous trees. Additionally, Richardson and O'Keefe (2009) show that the GC models are the best for 20 of 33 species, while the GDD only models are the best for 13 of 33 species at Harvard forest. Chiang and Brown (2007) also show that the phenology simulations from a three-parameter GC model suite the observations for all 17 selected species at Harvard forest. Thus, thorough tests of budburst models, including and excluding the chilling requirement, demonstrate that GC models work well for broadleaf deciduous temperate species.
Previous GC model comparisons (e.g. Richardson et al 2012) and sensitivity experiments (e.g. Jeong et al 2012) have mostly focused on deciduous plants in temperate regions, with less attention to temperature limited high latitude ecosystems. The Arctic has experienced severe changes, mainly characterized by an amplification of surface air temperature (Serreze et al 2009, Screen and Simmonds 2010, Screen et al 2013, which is altering the phenology of Arctic vegetation, favoring earlier leaf-out and photosynthesis (Randerson et al 1999, Euskirchen et al 2006, Zeng et al 2011, Piao et al 2013, Melaas et al 2018. This study explicitly investigates phenology modeling of Arctic deciduous types. It is suggested that the sensitivity of leaf budburst to accumulated growing degrees and chilling temperatures varies among different deciduous plant types in arctic tundra and boreal ecosystems (Euskirchen et al 2014, Delpierre et al 2016. In addition, post-snow melt air temperature is suggested as the dominant factor controlling leaf budburst and growth phenology in Arctic tundra and boreal ecosystems (Delbart andPicard 2007, Wipf andRixen 2010). Delpierre et al (2016) also demonstrate that for a majority of boreal tree species, a certain duration of chilling is needed for the release from endodormancy, followed by leaf budburst. Thus, we hypothesize that all the Arctic deciduous species need certain exposure to cold temperature for budburst, which regulates the vegetation phenology and carbon fluxes. We focus our analysis on Alaska, where aircraft and in situ carbon flux observations are available for the GC model benchmark. The specific objectives of the study are to (1) estimate the sensitivity of leaf budburst and photosynthetic onset of varied Alaskan deciduous vegetation types (e.g. deciduous forest, deciduous shrub, and C3 arctic grass) to a GC model and (2) where 'day' represents Julian day and T is surface air temperature. This scheme predicts that leaf budburst occurs when where NCD is defined as the number of chill days (any day with daily mean surface air temperature (T AIR ) less than 5°C constitutes a chill day) and T Jan-Feb is daily mean T AIR (°C) averaged over January and February (Jeong et al 2013). We define the right side of equation (4) as the GC threshold (figure 1). Here, a, b, c, and d are adjustable parameters (table S1 is available online at stacks.iop.org/ERL/15/025006/mmedia). We find that the ranges of a, b, and c developed in Jeong et al (2012) can well represent temperate deciduous forest spring phenology but some of the parameter groups cannot accurately represent the spring phenology in Alaska. Thus, we re-evaluate the adjustable parameters in Jeong et al (2012) and select two sets of coefficients (i.e. exp1 and exp2) that reasonably represent the spring phenology of Alaskan deciduous vegetation, as listed in table S1. Leaf growth phenology of tundra in North Slope Alaska is postponed more by exp2 than by exp1 (figure not shown; text S1). Thus, we design the third experiment by applying the parameters of exp1 to broadleaf deciduous shrub, which is the dominant vegetation type over North Slope Alaska as suggested by CLM plant functional type (PFT) maps (section 2.1.3), and applying the parameters of exp2 to broadleaf deciduous boreal tree type and C3 arctic grass (figure S1). This experiment is denoted exp3 (table S1).

The surface vegetation map of CLM
The PFT map in CLM is developed by using moderate resolution imaging spectroradiometer (MODIS) land surface mapping (Lawrence and Chase 2007). In Alaska, the vegetation cover fraction for needleleaf evergreen tree, broadleaf deciduous tree, broadleaf deciduous shrub, and C3 arctic grass is 21.2%, 3.6%, 36.5%, and 19.1%, respectively (table S2). We acknowledge that mapping the vegetation communities over the pan-Arctic and applying them to TBMs remain a challenge. Compared to the current CLM PFT map, a more diverse vegetation functional classification, such as a species-group based classification (Euskirchen et al 2014), is needed but beyond the scope of this study. Further, vegetation mapping over tundra, such as the North Slope Alaska, is more uncertain (Berner et al 2018). Nevertheless, the CLM PFT map prescribes shrubland over North Slope Alaska representing low-lying shrubs, consistent with satellite-observation and regression-algorithm based shrub maps (Beck et al 2011, Berner et al 2018, in which low-lying wetland areas on the coastal plain of the North Slope have the lowest shrub cover and ∼50% of the North Slope areas have shrub cover fraction 50% (figure S1).

Model setups
We use the atmospheric forcing data of Climatic Research Unit-National Centers for Environmental Prediction Version 7 (CRUNCEP7) at 0.5°× 0.5°s patial resolution as the forcing data of CLM (Lawrence et al 2019). Our primary analysis, which focuses on the sensitivity of leaf budburst and photosynthetic onset to GDD and chilling requirement, uses the default CLM and CLM with the GC model integrated (i.e. CLM_GC) driven by CRUNCEP7 and the surface PFT map from CESM1.2. Here, CLM and CLM_GC are run at the 0.9°latitude ×1.25°longitude and 30 min spatiotemporal resolution during 2007-2016, Figure 1. The diagram of GDD (the red line; equation (3)) and GDD plus chilling requirement (the navy line; GC threshold defined by the right side of equation (4)) regulated leaf budburst and GPP onset. The leaf budburst starts when GDD is larger than the GC threshold, and the red line and blue line cross at the budburst starting point. The GPP and surface air temperature (T AIR ) data are obtained from site measurements at Imnavait Creek watershed (68.4°N, 149.2°W; Euskirchen et al 2016). GDD and GC threshold in equation (4) are calculated with the T AIR data at this site. In this example, there is not a day with daily T AIR higher than 5°C when observed leaf budburst starts.
when both CRUNCEP7 (1901-2016), a remotesensing data-driven product (text S2), and the in situ observational carbon fluxes (text S3) are available. We use the spun-up surface conditions (i.e. with spun-up surface carbon and nitrogen pools) provided by NCAR to initialize the model. The model runs are summarized in table 1. (3) environmental effects using air temperature, soil temperature, and downward shortwave radiation from the North American Regional Reanalysis (NARR; Mesinger et al 2006). This PVPRM based and further constrained GPP product (herein PVPRM-SIF GPP) is validated against eddy covariance data and shows high consistency (Parazoo et al 2018). The constraint (1) above is also used to constrain ecosystem respiration (ER) from PVPRM, and the constrained ER is used together with the PVPRM-SIF GPP to obtain NEE. This NEE is further optimized using the Weather Research and Forecasting-Stochastic Time-Inverted Lagrangian Transport framework (Lin et al 2003) and atmospheric CO 2 vertical profiles obtained in the lower atmosphere across Alaska from April to November during Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE) campaigns ( The validations of PVPRM-SIF GPP and CARVE-OPT NEE provide confidence in our ability to benchmark CLM carbon dynamics in deciduous ecosystems. We evaluate CLM regional simulations during 2012-2014, and more detailed information on these two datasets is in text S2.

Surface vegetation classes over Alaska
PVPRM surface vegetation distribution are determined by the tree, shrub, and tundra classes over Alaska with the mean cover fractions 25%, 16%, and 44%, respectively; deciduous trees represent a very small percentage (1%) of Alaskan vegetation cover (table S4). PVPRM uses seven vegetation classes over Alaska (table S4), while CLM uses five PFTs (table S2). The primary difference between PVPRM class and CLM PFT maps is that the dominant vegetation in North Slope Alaska (68°N and north) is described as tundra by PVPRM and represented as a combination of shrub and C3 arctic grass in CLM PFT (figures 2(a) and S1).
To evaluate if CLM can reasonably simulate the vegetation phenology over Alaska, we carry our analysis over three vegetation groups by using the vegetation class map from PVPRM (figure 2). According to the bulk freeze/thaw status of the Alaskan land surface (text S4), the landscape thawing day of North Slope Alaska is typically ∼10-30 d later than that of southern Alaska (figure S3). Thus, we classify the vegetation class in North Slope Alaska as Alaskan northern tundra, separating it from southern Alaskan tundra (figure 2(a)). We combine the southern regions covered by either shrub or by shrub tundra (table S4) into Alaskan shrub ( figure 2(b)). According to both PVPRM and CLM, deciduous boreal forest represents a small fraction of Alaskan land cover. Considering the large model grid size (0.9°latitude ×1.25°longitude), deciduous boreal forest cannot be studied separately from CLM grid cells. Thus, we treat the three forest classes (table S4) as a single vegetation group, named as Alaskan forest (figure 2(c)). We identify one dominant vegetation group in each PVPRM and CLM grid cell. Vegetation is considered dominant when one of the three vegetation cover fractions exceeds 50% within a grid cell. Overall, our analyses are carried out separately over Alaskan northern tundra, Alaskan shrub, and Alaskan forest.

The definition of budburst start and GPP onset
In the default CLM, carbon stored in storage pools will be transferred to the display pools (e.g. leaf carbon pool) when GDD exceeds the GDD threshold (i.e. phenology onset). Leaf carbon is a function of LAI, and the phenology onset triggers the start of LAI accumulation (Oleson et al 2013). At this stage, the model simulates small GPP values and slow GPP growth in the Arctic, since the low temperature slows the vegetation growth. Thus, we define the leaf budburst date as the mean day of year (DOY) when GPP is between 1% and 10% of the peak value of annual GPP (i.e. GPP MAX ; figure 1). We also define the GPP onset date as the mean DOY when GPP is between 10% and 20% of GPP MAX for that year (Parazoo et al 2018). This definition can account for observation noise and range of transition dates from slow to rapid spring recovery in arctic tundra and boreal ecosystems. The number of days between GPP onset and GPP MAX is defined as growth season length. Both the definitions of leaf budburst and GPP onset are tested by using site observations at Imnavait Creek watershed (68.4°N, 149.2°W; Euskirchen et al 2016; text S5).

PVPRM-SIF inferred Alaskan GPP onset
We investigate how the factors in equation (4) affect dates of leaf budburst for each vegetation group. We first define a 'latest day of budburst', denoted as D BUDBURST, LAST (table 2), by searching for latest leaf budburst dates in PVPRM-SIF over the period 2012-2014 for each vegetation group. We then use T AIR from NCEP NARR, the meteorological forcing used to obtain PVPRM-SIF GPP, and the parameter values of exp1 to calculate NCD (number of chill days), T Jan-Feb (the mean T AIR in January and February), GC threshold (values of the right side of equation (4)), and GDDs for all days from 1 January to D BUDBURST, LAST in each Alaskan vegetation group (table S5). These calculations are used to investigate how the factors in equation (4) affect PVPRM-SIF GPP suggested DOY of leaf budburst in the three Alaskan vegetation groups.

CLM simulated Alaskan GPP onset
We first investigate if leaf budburst is regulated by T Jan-Feb and NCD by using PVPRM-SIF GPP and NARR T AIR (equation (4)). Based on section 2.3.4, we find that GDD and chilling exposure can explain the timing of leaf budburst in Alaskan northern tundra and Alaskan shrub (text S6). Since T AIR is the primary factor regulating leaf budburst and GPP onset in highlatitude regions (Euskirchen et al 2014), we compare  the T AIR anomalies from CRUNCEP7 with that from NARR during March, April, May, and June over Alaskan northern tundra, shrub, and forest to justify our use of CRUNCEP7. The results suggest that the interannual variability of CRUNCEP7 T AIR matches that of NARR T AIR , and it is suitable to study the variability of leaf budburst and GPP onset in CLM driven by CRUNCEP7 (text S7 and figures S4 and S5). We then use PVPRM-SIF to benchmark CLM estimated GPP onset during 2012-2014. CLM suggested mean DOY of GPP onset are 3, 4, and 1 d earlier than that suggested by PVPRM-SIF in Alaskan northern tundra, shrub, and forest, respectively. Compared to PVPRM-SIF, CLM suggests early DOY of GPP onset ranges from 1 to 12 d across the three vegetation groups (table 2). The spatial difference of GPP onset between CLM and PVPRM-SIF (CLM minus PVPRM-SIF) also shows that CLM suggests ∼38% of Alaska has earlier DOY of GPP onset (figure S6). By calculating the root mean square error (RMSE) of the differences between PVPRM-SIF and CLM simulated GPP, we find that CLM_GCexp1 simulated GPP onset has the least divergence (RMSE=1.7) from that of PVPRM-SIF among all the CLM simulations (table 2). Thus, CLM_GCexp1 gives a best estimation of DOY of GPP onset in the three ecosystems. CLM_GCexp1 postpones the DOY of GPP onset by 5±8, 6±8, 1±6 d compared to CLM and reduces the areas with early GPP onset to ∼29% of Alaska.
We further study the impacts of the GC model on DOY of GPP onset in the three Alaskan vegetation groups. In Alaskan northern tundra, CLM GPP onset is 2-9 d earlier than in PVPRM-SIF GPP, with the largest difference (i.e. 9 d) in 2014 (table 2). By integrating the GC model with different parameter groups (table S1) into CLM, GPP onset is postponed by a range of 2-20 d over the span of experiments and years. Even though 2014 has the warmest March, April and May (MAM) in Alaskan northern tundra, the delayed GPP onsets are best simulated by CLM_GCexp1 (i.e. 2 d later than the GPP onset of PVPRM-SIF) in 2014. Among the three parameter experiments, CLM_GCexp1 has the best representation of GPP onset compared to PVPRM-SIF, with delays ranging from 1 to 4 d (table 2). In Alaskan shrub, CLM GPP onset is 1 d earlier in 2012, 1 d later in 2013, and 12 days earlier in 2014 than in PVPRM-SIF. This result suggests that without including the GC model, CLM GPP onset diverges from that of PVPRM-SIF during the warm MAM of 2014. Compared to CLM, CLM_GCexp1 gives an improved estimation of GPP onset in Alaskan shrub; in particular, GPP onset is postponed by 13 d in 2014, which is much closer to PVPRM-SIF (table 2). Thus, the GC model also improved the GPP phenology over Alaskan shrub, especially in the years with warmer MAM T AIR (figures S4(b) and (e)). The GPP onset difference between PVPRM-SIF and CLM simulation in Alaskan forest is less than 2 d during 2012-2014 (table 2). In CLM, needleleaf evergreen boreal trees uses an evergreen phenology model, in which the leaf phenology depends on leaf longevity. Due to the large cover fraction of needleleaf evergreen boreal trees in the surface PFT map of CLM (figures S1(a) and (b), table S2), CLM_GC does not show much sensitivity of GPP onset to the chilling requirement. CLM_GCexp1 can exactly represent GPP onset dates in these three years over Alaskan forest. We also include site-level comparisons based on the measurements of three tundra types at one tundra site and one forest site in Alaska. The comparisons show that CLM simulates earlier GPP onset at both the tundra and forest sites, and the GC model postpones the GPP onset at arctic tundra sites (text S3).
The GC model is also crucial to the growth season length of Alaskan vegetated regions. During 2007-2016, CLM simulated mean DOY of GPP onset is 150±4, 138±4, and 134±5 in Alaskan northern tundra, shrub, and forest, respectively, and the growth season length for the same three ecosystems is 41±6, 32±10, and 34±7 d. CLM_GCexp1 postpones the mean DOY of GPP onset to 155±3, 142±4, and 135±4, and shortens the growth season length by 12%, 13%, and 4% in these three ecosystems ( figure 3). During 2008-2016, the DOY of GPP onset at the Imnavait Creek watershed (text S3) is 156±5, 146±3, and 158±5 as suggested by the site measurements, CLM, and CLM_GCexp1, respectively. CLM_GCexp1 reduces CLM simulated growth season length (35 d) by 35%. The DOY of leaf budburst has the similar variability (i.e. standard deviation) values to the DOY of GPP onset (figure not shown). Further, these variability numbers of CLM and CLM_GCexp1 show that CLM_GCexp1 only slightly alters GPP onset variability, and the DOY of GPP onset map affirmed this conclusion ( figure 4). Additionally, the spatial variability of DOY of GPP onset indicates that Alaskan northern tundra has smaller variability than other Alaskan vegetated regions. Since CRUNCEP7 T AIR always shows a larger variability in Alaskan northern tundra than in other Alaskan regions in April, May and June (AMJ; figure not shown), we infer that Alaskan shrub and forest (mixed with broadleaf deciduous tree, broadleaf deciduous shrub, and C3 arctic grass) have a larger sensitivity to T AIR variations than Alaskan northern tundra (dominated by broadleaf deciduous shrub) as suggested by CLM. The physiological reason for this difference is beyond the scope of this study.

The GC model induced carbon flux changes in Alaskan spring
Besides GPP phenology, the AMJ carbon fluxes are also affected by the GC model (Jeong et al 2012). Thus, we quantify the impacts of the GC model on AMJ carbon flux simulations. According to section 3.1, CLM_GCexp1 has a best estimation of GPP onset in all the three Alaskan vegetation groups with respect to PVPRM-SIF. Thus, we further study the impacts of CLM_GCexp1 on AMJ GPP and net ecosystem production (NEP; negative value of NEE) during 2012-2014.
In Alaskan northern tundra, AMJ GPP is 0.4 and 0.68 g C m −2 d −1 for PVPRM-SIF and CLM, respectively. By carrying out CLM_GCexp1, AMJ GPP decreases by 0.02 g C m −2 d −1 (table 3). Likewise, in  Alaskan shrub, AMJ GPP is 0.84 and 1.89 C m −2 d −1 for PVPRM-SIF and CLM, respectively, and CLM_GCexp1 leads to a GPP reduction of 0.39 g C m −2 d −1 (14%) over default CLM (table 3). We find a similar pattern of GPP reduction associated with CLM_GCexp1 over Alaskan forest (table 3). Here, CLM overestimates AMJ GPP in all three Alaskan vegetation groups compared to PVPRM-SIF GPP, especially over Alaskan shrub. CLM_GCexp1 simulated GPP fluxes suggest that among the three Alaskan vegetation groups, Alaskan shrub has the largest fraction of AMJ GPP reduction with the value of 14% (table 3). We also investigate AMJ GPP flux over all the Alaskan vegetated region (vegetation cover fraction 50%). The results also show that AMJ GPP is overestimated by CLM (by 110%) compared to PVPRM-SIF. CLM_GCexp1 reduces the AMJ GPP by 12% (0.28 g C m −2 d −1 ; figure 5(a)). Thus, the GC model could reduce the high AMJ GPP bias in CLM.
We use CARVE-OPT NEP to further benchmark CLM. Integrating CARVE-OPT NEP over all Alaskan vegetated areas, we find a mean AMJ uptake during 2012-2014 of 0.17 g C m −2 d −1 , indicating a growth season carbon sink. In comparison, the mean AMJ NEP in Alaskan northern tundra, shrub, and forest is −0.47 g C m −2 d −1 , 0.17 g C m −2 d −1 , and 0.61 g C m −2 d −1 , respectively, indicating the Alaskan growth season sink is associated with shrub and forest (table 3). Compared to CARVE-OPT, CLM overestimates net Alaskan AMJ uptake by a factor of two (NEP=0.38 g C m −2 d −1 ), producing sinks in all three vegetation groups (NEP=0.06, 0.41, and 0.59 g C m −2 d −1 in northern tundra, shrub, and forest, respectively). CLM_GCexp1 reduces this sink by 0.22 g C m −2 d −1 ( figure 5(b)), in closer agreement with CARVE-OPT. The variations of Alaskan carbon fluxes at the yearly time scale are discussed in Text S8.

Discussion
Model intercomparisons have repeatedly shown an early bias in predicted photosynthetic-growing season onset in northern high latitudes on local to regional scales (Richardson et al 2012, Peng et al 2015, Commane et al 2017. Previous studies have hinted that leaf budburst and subsequent GPP onset in temperature limited ecosystems are sensitive to accumulated GDD as well as cold temperature (chilling) exposure. By integrating a new phenology model into CLM that accounts for both effects (CLM_GCexp1), we are able to successfully postpone GPP onset in Alaskan northern tundra and shrub by 5±7 and 4±8 d and reduce spring GPP magnitude by 2.5 g C m −2 and 32.1 g C m −2 during 2007-2016. These variations lead to a reduction in spring (AMJ) GPP of 10% from 2007 to 2016. Overall, we find CLM_GCexp1 delays GPP onset by 2-13 d compared to the unadjusted CLM over Alaska from 2012 to 2014, in closer agreement with an observationally Table 3. PVPRM-SIF, CLM, CLM, and GCexp1 suggested 2012-2014 mean AMJ GPP (g C m −2 d −1 ) and NEP (g C m −2 d −1 ) values over Alaskan northern tundra, Alaskan shrub, and Alaskan forest.

PVPRM-SIF CLM CLM_CGexp1
Alaskan  constrained benchmark (PVPRM). The model simulated GPP in Alaskan forests is not as strongly affected in the new model due to dominance of needleleaf evergreen tree over Alaska. These results suggest that the thermal forcing models with a chilling requirement have the potential of reducing the uncertainty of TBM simulated spring phenology in high-latitude regions.
The timing of leaf budburst is mostly unchanged in CLM_GCexp1, with the exception that leaf budburst in tundra is further delayed relative to PVPRM-SIF from 1 d in CLM to 8 d in CLM_GCexp1 over 2012-2014 (table 2). We note a couple caveats here. First is the transition to spring green-up is very slow and noisy from PVPRM (and in situ observations), making detection of a ∼5% threshold for leaf budburst uncertain. Second is that PVPRM-SIF uncertainty increases moving into earlier spring with reduced signal and increased noise in the spaceborne SIF constraint. As such, further testing of budburst in Alaskan deciduous vegetation types is needed.
The new model also alters the structure of peak growing season GPP. Specifically, CLM_GCexp1 shifts the timing of peak GPP by −6, 7, 3 d and amplifies the magnitude of peak GPP by 4.7%, 0.9%, and 9.2% over Alaskan northern tundra, shrub, and forest, respectively, during 2007-2016 (figure 3; 3.6% increase on average). It is likely that the increase of annual GPP peak could be associated with a postponed leaf budburst, which might allow more nitrogen accumulation in soil during and after soil thaw and then trigger a higher GPP amount after GPP onset (Larsen et al 2007, text S8).
To improve the consistency of annual maximum LAI values between MODIS and CLM, CESM Version 2 (CESM2.0) updates its high-latitude PFT map, producing a 26% increase in C3 arctic grass fraction in Alaska (table S2). CLM simulations based on this updated PFT map (herein CLM_P2) produce larger disagreement in GPP phenology (raised early bias) and magnitude (increased high bias) over Alaskan northern tundra and Alaskan shrub compared to simulations with the original map (table S5). These results and the improved vegetation map used by CLM_P2 imply that the functionality rather than the surface misrepresentation of PFTs induces the overestimated GPP in CLM. Thus, uncertainties still exist in GPP simulations in the photosynthetic module of CLM and other TBMs. Future studies at smaller spatial scales (i.e. Alaskan tower sites) are needed to identify the reasons, which could be associated with photosynthetic related prognostic variables (e.g. Vcmax) and nutrient limitations (Rogers et al 2019).
Parameter availability limits the application of the GC model to different PFTs in arctic tundra and boreal ecosystems over Alaska. We test all the parameters obtained from temperate forest observations in Jeong et al (2012) and (2013) and select the parameter groups that can be reasonably applied into CLM over Alaska.
All the three experiments, using the same model structure (equations (3) and (4)) but different parameter groups (table S1), show the sensitivity of 'GDD plus chilling requirement' to different parameter groups. Here, we find that with the same number of NCD, exp2 provides higher GC threshold values than exp1 (figure S7). Thus, exp2 requires a higher number of GDD suggested in equation (3), which takes a longer accumulative time period, to trigger the leaf budburst. This structure of the GC model explains the later GPP onset represented by exp2 than by exp1 (table 3). We also try to use the least-square curve-fitting method and the observations at the Imnavait Creek watershed to optimize the parameters. However, owing to the limited amount of data, the result does not coverage. Thus, acquiring and synthesizing more budburst records over the Arctic (Euskirchen et al 2014) and using varied approaches for parameter optimization are needed for carbon cycle studies and for the potential improvements to GPP phenology representation in CMIP models (Anav et al 2013) in high-latitude regions.
The existing studies on climate change and plant phenology in Arctic tundra are very limited compared to similar studies in lower latitudes (Diepstraten et al 2018). Additionally, the range of chilling temperatures needed for ending dormancy has not yet been identified for different Arctic deciduous vegetation types (Dantec et al 2014). This limitation is associated with the complexity of understanding the biological processes and the high cost in maintaining long-term observations to understand how these biological processes are affected by the changing environmental factors in Arctic tundra (Diepstraten et al 2018). Bjorkman et al (2015) show that the timing of snowmelt is also a strong driver of flowering phenology, especially for early-flowering species over Arctic tundra, while the flowering phenology for late-flowering species depends on spring temperature. With the observational limitations, disentangling the impacts of different environmental factors on plant phenology and including these impacts into phenology representations of TBMs are still a challenge with respect to improving phenology simulations over Arctic ecosystems. Thus, high spatial and temporal resolution observations that can monitor vegetation phenology are needed to be advanced for carbon cycle studies in high-latitude regions.

Conclusions
We integrated a model which includes GDDs and a new chilling requirement into CLM. Compared to a model constrained by satellite observations (PVPRM-SIF), we find the chilling requirement can better represent photosynthetic onset in Alaskan deciduous vegetated regions. The revised model postpones photosynthetic onset over Alaskan northern tundra and shrub by 2-13 d over the period 2012-2014, with greatest improvement during warm springs (e.g. 2014). Our GC model also reduces the high AMJ GPP bias (compared PVPRM-SIF GPP) by ∼12% over vegetated Alaskan areas. Investigation of additional processes such as nitrogen limitation is needed to further reduce high GPP bias in CLM. This study represents a critical step forward in predicting Arctic deciduous vegetation phenology and its response to future warming.