Radiation‐constrained boundaries cause nonuniform responses of the carbon uptake phenology to climatic warming in the Northern Hemisphere

Abstract Climatic warming has lengthened the photosynthetically active season in recent decades, thus affecting the functioning and biogeochemistry of ecosystems, the global carbon cycle and climate. Temperature response of carbon uptake phenology varies spatially and temporally, even within species, and daily total intensity of radiation may play a role. We empirically modelled the thresholds of temperature and radiation under which daily carbon uptake is constrained in the temperate and cold regions of the Northern Hemisphere, which include temperate forests, boreal forests, alpine and tundra biomes. The two‐dimensionality of the temperature‐radiation constraint was reduced to one single variable, θ, which represents the angle in a polar coordinate system for the temperature‐radiation observations during the start and end of the growing season. We found that radiation will constrain the trend towards longer growing seasons with future warming but differently during the start and end of season and depending on the biome type and region. We revealed that radiation is a major factor limiting photosynthetic activity that constrains the phenology response to temperature during the end‐of‐season. In contrast, the start of the carbon uptake is overall highly sensitive to temperature but not constrained by radiation at the hemispheric scale. This study thus revealed that while at the end‐of‐season the phenology response to warming is constrained at the hemispheric scale, at the start‐of‐season the advance of spring onset may continue, even if it is at a slower pace.

, successional transition (Chuine, 2010) and plays a role in the feedback between vegetation and climate (Peñuelas & Filella, 2009;Richardson et al., 2013). The warming-induced lengthening of the growing season has increased the carbon uptake (Le Quéré et al., 2009), offsetting atmospheric carbon from human emissions. However, how phenology may respond to future warming and whether vegetation will increase carbon sequestration remains unclear , which adds uncertainty to future atmospheric carbon concentration and, thus, climate projections.
Winter chilling, forcing requirements, and photoperiod trigger the start of foliar phenology (Richardson et al., 2013). Previous research suggests that chilling requirements will constrain the advance of leaf unfolding in deciduous forests by reducing temperature sensitivity (Fu et al., 2015), whereas other factors such as precipitation (Peaucelle et al., 2019) and photoperiod (Körner & Basler, 2010;Meng et al., 2021;Zohner et al., 2016) may increase the heat requirements during the ecodormancy stage and, as a result, slow down the warming-induced advance of the leaf unfolding. With regards to leaf senescence, a recent study found that increased productivity during the growing season counteracts the warming-induced delay in leaf senescence (Zani et al., 2020). These findings are based on the study of leaf phenophases in deciduous forests obtained from in situ observations (Templ et al., 2018) or remotely-sensed vegetation greenness indices (Julien & Sobrino, 2009), such as the normalized vegetation difference index (NDVI).
The key conceptual framework in these studies is that conditions preceding the phenological dates affect the timing of tree phenophases (Chuine et al., 2013). Phenology modelling has largely accounted for chilling and forcing requirements in decisuous forests (Fu et al., 2020). However, the start and end of the photosynthetically active season are directly influenced by current meteorological conditions, not by conditions preceding the phenological dates. This is evidenced by a decoupling between remotely sensed vegetation greenness and proxies of photosynthetic activity (Jeong et al., 2017;Yin et al., 2020;Zhang, Commane, et al., 2020), which indicates that vegetation might present leaves but these are not photosynthetically active because meteorological conditions are restricting photosynthetic activity at the moment. For instance, vegetation productivity declines in agreement with the decrease in radiation intensity during autumn (Zhang, Commane, et al., 2020), and a lack of available light may, thus, prompt the end of the growing season. Similarly, water availability determines the start and end of the photosynthetically active season in tropical dryland ecosystems (Eamus & Prior, 2001;Jiao et al., 2021;Zhang, Parazoo, et al., 2020).
Given the link between photosynthesis and current meteorological conditions, the seasonality of carbon uptake can be modelled with meteorological variables, and the start and end of the growing season can be determined by constraint functions (Jolly et al., 2005).
The constraint functions define the photosynthesis-inhibiting thresholds that restrict the carbon uptake during the growing season. Using this phenology modelling framework, carbon uptake phenology can be, thus, determined by the most limiting factor (Zhang, Parazoo, et al., 2020). The start and end of the season occur when a limiting factor, such as air temperature and radiation intensity exceeds the threshold under which photosynthesis is constrained.
One limiting factor of the phenology of photosynthesis is the low light availability (Zhang, Commane, et al., 2020). Photosynthetic activity as a function of radiation, with constant temperature and soil water content conditions, is equivalent to light-response curves (Herrmann et al., 2020), which represent the restrictions of photosynthetic activity by low levels of radiation, even if temperature, soil water content, and other factors influencing photosynthesis are favourable. The restriction is due to the inoperability of chloroplasts and light-dependent reactions in the absence of photosynthetically active radiation (Johnson et al., 2008).
Here, we aimed to clarify the role of temperature and radiation constraints during the start and end of the photosynthetically active season and how these constraints change among ecosystems and over space. To achieve this, we proposed a phenology model based on the law of limiting factors that predicts the start and end of the photosynthetically growing season. We evaluated the current temperature and radiation constraints on the phenology of temperate and cold regions of the Northern Hemisphere, where temperature and radiation mostly regulated the vegetation phenology . Lastly, we estimated the potential lengthening of the growing season in a warming scenario in which radiation might take over the role as a limiting factor of photosynthetic activity phenology.

| Data
We used daily in situ records from 63 sites of the FLUXNET2015 Tier 1 data set (Pastorello et al., 2020) (see Figure S1 for the locations of the sites). The biome types at the FLUXNET sites were characterized using the MODIS MCD12Q1 V6 product (Friedl & Sulla-Menashe, 2015) for the mode of land cover during 2001-2019 and the RESOLVE Ecoregions 2017 map (Dinerstein et al., 2017). We selected the sites with at least 4 years of record and categorized as deciduous broadleaved forests (DBF; 12 sites), evergreen needleleaved forests (ENF; 25 sites), mixed forests (MX; 7 sites), and grasslands and wetlands (GRA and WET; 19 sites) in the tundra, boreal, and temperate biomes of the Northern Hemisphere according to the RESOLVE Ecoregions 2017 map. These FLUXNET sites represent 643 site-years and are located >30° N. The sites cover different range of years from 1991 to 2014. We used daily gap-filled GPP reference from GPP versions using model efficiency (MEF) and obtained from the day-time partitioning method (Pastorello et al., 2020).
For the climatic data, we used daily averages of air temperature and incoming shortwave radiation at 0.1 arc degrees from the ERA5-Land hourly data (Muñoz-Sabater et al., 2021). The phenology dates were extracted from the 4-day clear-sky daily contiguous solarinduced chlorophyll fluorescence (SIF) estimates at 0.05° from the OCO-2 (Orbiting Carbon Observatory-2) (Zhang et al., 2018) (Sun et al., 2017). The temperate and cold regions include tundra, boreal forests, temperate broadleaf and mixed forests, and temperate coniferous forests in the RESOLVE Ecoregions 2017 map ( Figure S1).

| Phenology extraction
We estimated the Start of Season (SoS) and End of Season (EoS) using the maximum-separation method . The maximum-separation method is a threshold-based method that can effectively estimate phenological metrics without the need of time series pre-processing, which tends to distort the seasonality of the time series. The algorithm first runs a moving window that calculates the proportion of observations that are above a given threshold th before the central day of the moving window (p before ) and after the central day (p after ). A moving window greater than 90 days and lower than the length of the growing season is advised for the method to be effective . In this study, the length of the moving window was 120 days. Then, the algorithm computes the difference between the proportions (p before -p after ) for each day of year in the time series. The SoS is calculated as the date when the difference in proportions (p before -p after ) reaches a minimum, and the EoS is calculated as the date when the difference is maximal. The The I max and I min are the maximum and minimum values of the time series (e.g. FLUXNET GPP or the growing-season index). The parameter P ranges from 0 to 1. We set P equal to 0.2, which represents a percentage of 20% of the amplitude. This low percentage ensures the reliable detection of the first and last stages of the SoS and EoS. The maximum-separation method can be applied to binary time series, which is the case of the climate-constraint model (See Section 2.4. Climate-constraint model). In a binary time series in which values are 0 or 1, the SoS and EoS would be the same for any percentage P greater than 0 and lower than 1. For this reason, we used a constant threshold th equal to 0.5 for the climate-constraint model.

| Phenology modelling
We proposed a phenology model that incorporates the constraints of temperature and radiation on the vegetation productivity. The model first generates a binary time series where daily observations are categorized as 0 (vegetation productivity is inhibited or constrained) and 1 (vegetation productivity is active or unconstrained by temperature and radiation). In the second step of the model, the start and end of season dates are extracted from the binary time series using the maximum-separation method. We named the model as climate-constraint model because climate time series are used to generate a binary time series of constrained and unconstrained daily observations. The binary time series is analogous to the growingseason index (Jolly et al., 2005), which is generated from cut-off functions from daily temperature, daylength, and vapor-pressure deficit (VPD), and ranges from 0 and 1 (0 = vegetation is constrained and 1 = vegetation is unconstrained).
The modelled SoS and EoS were compared with estimated SoS and EoS extracted from in situ GPP time series at the FLUXNET sites. The performance of the model was evaluated in terms of root-mean-squared-error (RMSE). As the growing-season index, the climate-constrained model assumes no forcing requirements for the SoS and EoS; the growing season starts as soon as climate conditions are favorable for vegetation growth, respectively.
To test this assumption, we compared the performance of the climate-constrained model with the photo-threshold model (Meng et al., 2021), a model based on the growing-degree-day model. We also compared the performance of the SoS and EoS extracted from the growing-season index.

| Climate-constraint model
The approach for phenological modelling categorizes daily time series of temperature T(t) (°C) and incoming shortwave radiation R(t) (W m −2 ) into a binary time series B(t) with values 0 (constrained) and 1 (unconstrained). The classification was done with a threshold function (Equation 2); vegetation was considered unconstrained when daily temperature and radiation were above the threshold function. The threshold function has the form of a rational function with one asymptote for temperature (T thresh ) and another for radiation (R thresh ) ( Figure 1a). The curvature of the rational function is defined by parameter C (unitless). The parameter C in Equation (2) reflects the covariation of temperature and radiation constraints. Instead of a constant threshold for temperature and radiation, the parameter C reflects the curvature of the threshold function and provides a threshold that depends on the levels of temperature and light availability (i.e., high levels of temperature reduce the radiation constraint). We used a rational function owing to the relationship of FLUXNET GPP with air temperature and shortwave radiation ( Figure 1b).
The time series of FLUXNET GPP and ERA5-Land air temperature and shortwave radiation were used to calibrate the parameters T thresh , R thresh , and C. To accomplish this, the FLUXNET GPP was converted to a binary time series: 0 when GPP was < GPP thresh and 1 when GPP was > GPP thresh C m −2 d −1 . We used a GPP thresh equal to 2 C m −2 d −1 , which is approximately the median GPP of the 192,709 observations and represents a low GPP threshold that ensure the first and last stages of the SoS and EoS, respectively. The calibration was done using the simplex search method (Lagarias et al., 1998). EoS by the maximum-separation method.

| Growing-season index
The growing-season index was designed to reproduce the NDVI time series from air temperature, daylength, and VPD (Jolly et al., 2005).
The growing-season index is calculated from cut-off functions for F I G U R E 1 Graphical representation of the calculation of parameters in the climate-constraint model. (a) Graphical representation of the threshold function (black line) and the angle θ. The threshold function consists of a rational function defined by the asymptotes in T thresh and R thersh . The angle θ pheno represents the radiation constraint of temperature and radiation during a phenological event. (b) Mean gross primary production (GPP) as a function of air temperature and incoming shortwave radiation for all observations of the FLUXNET sites in the cold and temperate regions of the Northern Hemisphere.
parameters of the cut-off functions. However, we calibrated these parameters for our case study because (1) we used different climate data, (2) we used mean air temperature instead of minimum air temperature, and (3) we aimed to reproduce the FLUXNET GPP, representing vegetation productivity, instead of the NDVI, which is an indicator of green biomass. The calibration of the parameters was done using the simplex search method.
The growing-season index can be used as a phenology model. As the climate-constrained model, phenological dates can be extracted from the growing-season index. We used the maximum-separation method and a dynamic threshold of 20% of the amplitude to determine the SoS and EoS dates.

| Photo-threshold model
The photo-threshold works similarly as the growing-degree-day model (Hänninen, 2016) but replaces the arbitrary date when forcing starts accumulating (generally January 1) with the date when daylength reaches the threshold DL start (Meng et al., 2021).

| Estimation of temperature and radiation constraints on phenological dates
To quantify the limitations of radiation and temperature on the start and end of the growing season, we reduced the two-dimensionality of the temperature-radiation constraint to one single variable; θ.
The angle θ was calculated geometrically (Equation 9) based on the threshold function (Equation 2) and given the points Pa, Pb, and Pc, where Pa and Pb have fixed coordinates in (T thresh , R thresh ) and (25, R thresh ), respectively. The coordinates of Pc are given by the temperature and shortwave radiation at the time of the phenological event (T pheno and R pheno ). We calculated T pheno and R pheno as the mean temperature and radiation 7 days before and after the phenological event (SoS and EoS) to avoid noisy values of θ due to variation on daily temperature and radiation. A graphical representation of the calculation of θ is shown in Figure 1a. The temperature values (T pheno and T thresh ) were scaled by dividing by 25°C, and radiation values (R pheno and R thresh ) were scaled by dividing by 300 W m −2 . These scaling values represent approximately the 95th percentile of all temperature and radiation observations in the FLUXNET sites.
A high temperature and low amount of radiation at the time of the phenological event lead to low values of θ proximal to 0°, and a low temperature and high amount of radiation at the time of the phenological event lead to high values of θ proximal to 90°. Values of θ <0 or >90° were clamped to 0 and 90°, respectively.

| Estimation of the potential lengthening of the growing season
We converted the current climate constraints on phenology into potential lengthening of the growing season. The potential SoS and EoS (SoS pot and EoS pot ) represent the dates under a warming scenario in which SoS and EoS would be restricted only by radiation.
This gives a practical quantification of the climate constraints on the phenological dates. For instance, a small difference between EoS and EoS pot indicates that future warming will not delay the EoS. For the calculation of SoS pot and EoS pot , we changed the threshold function in Equation 2 to a constant radiation threshold equal to R thresh .
The observation was considered inhibited when shortwave radiation was <R thresh and active when it was >R thresh . Analytically, this would correspond to the climate-constraint function (Equation 2) having a high temperature T(t) that is nearly infinity.

| Estimation of temperature dependency
The dependency of SoS and EoS on temperature was estimated for 2001-2020 using ERA5-Land and the phenological dates extracted from the SIF time series. The SoS and EoS were extracted using a dynamic threshold of 20% the amplitude. The mean air temperature during the phenological dates were calculated as the average daily air temperature during ±7 days the SoS and EoS. Temperature dependency was estimated using the coefficient of correlation between temperature and phenological dates obtained from SIF time series, which provides more spatial and temporal coverage (2001-2020) than the FLUXNET GPP records. We used Pearson correlations, which is termed temperature dependency, instead of the conventional metric of temperature sensitivity (Fu et al., 2015) because the latter overestimates the sensitivity when the temperature time series is highly variable (Keenan et al., 2020), which could potentially weaken the spatial analysis in our study. Pixels were masked if minimum winter radiation was below R thresh ; latitudes where winter radiation is above R thresh are potentially not limited by radiation and present SIF values greater than 0. The minimum winter radiation was estimated as the average shortwave radiation from day of year 345 to 365 and aggregated over the 2001-2020 period.

| RE SULTS
The parameters T thresh , R thresh , and C, which define the threshold function of temperature and radiation, were 2.56°C, 26.5 W m −2 , and 342 when considering all FLUXNET site-years (

| DISCUSS ION
Our results indicate that the photosynthetically growing season Our results highlight that temperature and radiation constraints differ spatially during the SoS and EoS, indicating that vegetation will respond to climate change differently depending on the region of the Northern Hemisphere. Radiation is already constraining the SoS and EoS in regions such as Europe, while the EoS is largely constrained by radiation in the Northern Hemisphere, as indicated previously (Zhang, Commane, et al., 2020), but not in highly elevated areas where temperature mainly limit photosynthesis in autumn.
Moreover, the advance of spring onset was overall not limited by radiation, and future warming can potentially lead to an advance of the spring onset of carbon uptake in evergreen needleaved forests, which respond promptly to favorable temperature and radiation conditions. F I G U R E 2 Comparison between observed and simulated start of the growing season (SoS) (in red) and end of the growing season (EoS) (in blue) at the FLUXNET towers in evergreen needleleaved forests (ENF), deciduous broadleaved forests (DBF), mixed forests (MX), and grasslands (GRA). Observed SoS and EoS were obtained from GPP time series from the FLUXNET towers. Simulated SoS and EoS were obtained from ERA5-Land temperature and radiation time series using the climateconstrained model. The titles of the charts report the mean error (ME), the root-mean-squared-error (RMSE), and the number of observations (n).

TA B L E 2
Root-mean-squared error (RMSE) between modeled and estimated phenological dates at the FLUXNET sites The climate-constraint model is similar to the growing-season index, which considers three separate cut-off functions for temperature, daylength, and VPD (Jolly et al., 2005). A difference between these models is that the climate-constraint model considers the interaction between the constraints of radiation and temperature, which might explain the better estimates of phenological dates. Our model used radiation instead of daylength. Daylength has been proposed as a factor controlling the spring onset (Körner & Basler, 2010) and leaf senescence (Way & Montgomery, 2015).

Climateconstraint model
Daylength also plays a role in regulating the forcing requirements in deciduous forests (Meng et al., 2021). However, radiation is more suited for modelling the start and end of the photosynthetically active season (Zhang, Commane, et al., 2020), primarily because photosynthesis is directly linked to light availability (Herrmann et al., 2020).
The cold-avoidance strategy of deciduous forests may explain their exception to the law of the limiting factors in the phenology of the carbon uptake, as revelead by the long difference between modelled SoS and flux-based SoS in deciduous forests. Deciduous trees respond differently than evergreen conifers to cold winters.
Most conifers in temperate and cold regions have a cold-tolerant strategy. Leaves in conifers resist cold, reject excessive radiation when temperatures are still low in the spring, and resume photosynthetic activity as soon as climatic conditions are favorable . In contrast, deciduous trees shed their leaves in autumn and leaf budbreak in spring (Chuine et al., 2013). This would account for our findings that deciduous trees remain dormant during spring, even when photosynthesis conditions are favorable. Such strategy would also avoid leafing out prematurely to prevent frost damage F I G U R E 3 Potential lengthening of the growing season in temperate and cold regions of the Northern Hemisphere.
(a) Magnitudes of the constraints on the start and end of season (SoS and EoS) by radiation and temperature. The temperature constraints increase with θ, and the radiation constraint decreases with θ. (b) Potential advance of the SoS (SoS − SoS pot ) and potential delay of the EoS (EoS − EoS pot ). SoS pot and EoS pot represent the start and end of the growing season if radiation was the only limiting factor. SoS and EoS were estimated using the OCO-2 SIF data set and temperature and radiation were extracted from the ERA5-Land for 2001-2020.

F I G U R E 4
Histograms of the coefficients (R) for the correlation between temperature (ERA5-Land) and phenological dates (estimated from the OCO-2 SIF time series) during 2001-2020. The asterisk represents that mean R was significantly different in the regions with the highest radiation constrain during the SoS (θ < 65) and in the regions with the highest temperature constrain during the EoS (θ > 50). The threshold of θ for the SoS represents the 5th percentile of all mean θ values and the threshold of θ for the EoS is the 95th percentile of the mean θ values for the 2001-2020 period. (Meng et al., 2021). The ecodormancy stage is only broken when a certain amount of heat is accumulated (Richardson et al., 2013).
Temperature was the main limiting factor during spring onset, but temperature and radiation have been limiting during autumn in temperate and cold regions of the Northern Hemisphere (Zhang, Commane, et al., 2020). observations have been made extensively over the last two decades (Pastorello et al., 2020;Zhang et al., 2018).
The asymmetric temperature variability over the season explains the different radiation and temperature constraints during the SoS and EoS. In turn, the different radiation and temperature constraints explain the nonuniform seasonal temperature dependency during F I G U R E 5 Potential lengthening of the growing season in the FLUXNET sites. Y-axis show the latitudinal gradients of the observed start of season (SoS; red dots) and the End of Season (EoS; blue dots) obtained from the GPP time series at the FLUXNET towers. The mean potential start of season SoS pot and end of season EoS pot (black dots) represent the Day of Year (DoY) when shortwave radiation surpasses the 27 W m −2 threshold. The DoY was offset such that the summer solstice (DoY solstice ) is the origin of the axis.
the SoS and EoS and across the temperate and cold regions, suggesting divergent responses of vegetation phenology to future climatic warming. The dependency of the SoS and EoS to temperature was significantly lower in regions where radiation constraints were the highest. These results indicate that the sensitivity of vegetation phenology to temperature decreases with the increasing constraints of radiation. Our results show that the senescence stage has a low temperature dependency due to the constraints of radiation, and this might be a reason for the lower magnitude in the EoS delay than the SoS advance (Menzel et al., 2006). In deciduous forests, other alternative mechanisms for leaf senescence include the growing-season productivity (Zani et al., 2020). SoS could continue its advance with future warming, even if it is at a lower pace (Fu et al., 2015), although the advance of spring onset would be highly spatially variable depending on the local radiation limitations. In needleaved forests, radiation should be taken into account when modeling how the phenology of the photosynthetically active season will change as the climate warms.

| CON CLUS IONS
This study showed that limiting factors determined the start and end of the photosynthetically active season in evergreen needleaved forests in temperate and cold regions of the Northern Hemisphere.
Our findings showed that temperature and radiation constraints on vegetation productivity differed depending on the regions of the Northern Hemisphere. In addition, the radiation constraints were different during the SoS and EoS, indicating that carbon uptake phenology will respond to future warming differently in spring than autumn. Temperature was overall limiting photosynthesis during the spring onset, but temperature and radiation were both constraining at the end of the growing season. The radiation constraint increases in autumn because temperatures remain high at that stage of the season, suggesting that future warming will not delay the end of the photosynthetically active season. In contrast, radiation limitation was low during the spring onset and the start of the season was highly sensitive to temperature, which indicates that the photosynthetically active season could continue overall its advance with future warming. Our results highlight the importance of taking radiation, temperature, and their covariance into account when modeling the photosynthetically active season in evergreen needleleaved forests and its responses to climatic warming.