# Decreasing rainfall frequency contributes to earlier leaf onset in northern ecosystems

Jian Wang, Desheng Liu, Philippe Ciais, Josep Peñuelas

► To cite this version:

Jian Wang, Desheng Liu, Philippe Ciais, Josep Peñuelas. Decreasing rainfall frequency contributes to earlier leaf onset in northern ecosystems. Nature Climate Change, 2022, 12, pp.386-392. 10.1038/s41558-022-01285-w. http://hal-03587397

HAL Id: hal-03587397

https://hal.science/hal-03587397v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Decreasing rainfall frequency contributes to earlier leaf onset in northern ecosystems

Jian Wang<sup>1,∗</sup>, Desheng Liu<sup>1,†</sup>, Philippe Ciais<sup>2</sup> and Josep Peñuelas<sup>3,4</sup>

Climate change substantially advances the leaf onset date (LOD) and regulates carbon uptake by plants. Unlike temperature, the effect of precipitation remains largely elusive. Here we use carbon-flux measurements, in situ records of leaf unfolding and satellite greenness observations to examine the role of precipitation frequency ($P_{freq}$, number of rainy days) in controlling the LOD in northern ecosystems (>30°N). Widespread decreases in $P_{freq}$ during the past three decades positively contributed to the advance in LOD, possibly due to increased exposure to radiation, exhibiting a dominant control of LOD over ~10% of the area. Lower $P_{freq}$ may also enhance chilling at night and warming at daytime, consequently leading to earlier LOD. We further develop a weighted precipitation growing-degree-day algorithm that projected a generally earlier LOD than currently predicted. These results highlight the need for a comprehensive understanding of the impacts of precipitation on LOD, which is necessary for improved projections.

The earlier leaf onset date (LOD) of northern vegetation under recent warming has been widely reported on the basis of eddy-covariance flux measurements<sup>1,2</sup>, in situ records<sup>3–6</sup> and satellite observations<sup>7,8</sup>. This shift in LOD can contribute to enhanced ecosystem productivity, with an earlier start of carbon uptake by plants<sup>1,9,10</sup>. Previous studies have focused mainly on the warming effect on LOD<sup>5,6,8</sup>, particularly in northern areas with a large carbon sequestration<sup>11,12</sup>. The impacts of precipitation on LOD, however, are largely elusive, partially because studies have focused on total amount of precipitation ($P_{total}$) without accounting for the frequency of precipitation ($P_{freq}$, number of rainy days)<sup>13,14</sup>. Exploring the impacts of $P_{freq}$ may therefore help us better understand the responses of LOD to climate change and reduce the considerable uncertainty in predicting LOD.

Recent warming has generally advanced spring LOD with a heterogeneous sensitivity to temperature (d°C⁻<fnref n="1" />) in northern ecosystems<sup>5,8</sup>. This is because the chilling accumulation (the amount of chilling received by plants during the first dormant stage—endodormancy) and heat requirement (the accumulated forcing temperature required for breaking the second dormant stage—ecodormancy) for budburst and leaf formation are controlled by temperature, precipitation, radiation and other forcings<sup>6,8,15</sup>. For example, it has been reported that an increase of daytime temperature by 1°C advanced satellite-based LOD by 4.7 days in Europe, 4.3 days in the United States and >10 days in northern Siberia and northwestern Canada during 1982–2011<sup>8</sup>. Unlike temperature, the effect of precipitation on LOD has received less attention due to complex mechanisms related to interactions with temperature, radiation, soil moisture and snow cover<sup>14,16,17</sup>. To date, $P_{total}$ has been used as the main characteristic of rainfall to look for influences on ecological processes and energy and carbon fluxes at terrestrial surfaces<sup>17–19</sup>. Extant studies suggested that an increase in $P_{total}$ may delay LOD in northern ecosystems<sup>14–16</sup> due to the increase in snowmelt heat requirement and the decrease in absorbed solar radiation. For example, larger winter precipitation acts as a critical cause of longer-lasting snow cover in high latitudes, leading to (1) lower temperature because of

increased snow-melting latent heat consumption and (2) a decrease of absorbed radiation due to high albedo of snow-covered surfaces<sup>15,16</sup>. Consequently, a wet winter could delay the heat accumulation required for leaf onset. Apart from $P_{total}$, $P_{freq}$ is crucial to assess climate change impacts<sup>20</sup>. On the basis of observations<sup>21</sup> and model projections<sup>22,23</sup>, $P_{freq}$ has been reported to be decreasing due to surface warming (thermodynamic contribution) and weakening of tropical circulation (dynamic contribution)<sup>24</sup>. Changes in $P_{freq}$ have notably affected plant growth and productivity by regulating run-off<sup>25</sup>, soil moisture<sup>26</sup>, exposure to high radiation and temperature, and energy fluxes<sup>27</sup>. Thus, interannual variations of $P_{freq}$ are expected to increase the effects on plant phenological transitions under warming, especially in arid regions. We hypothesize that changes in $P_{freq}$ control the effects of precipitation on LOD related to incoming radiation, heat and chilling accumulation and soil water availability. We tested this hypothesis by analysing gridded meteorological data, including near-ground mean temperature (T<sub>mean</sub>, °C), total cloudiness (C<sub>total</sub>), %, a proxy of solar radiation), $P_{total}$ (mm) and $P_{freq}$ (days), together with LOD proxies from four independent datasets at northern middle and high latitudes (>30°N): (1) 745 site-year records of gross primary productivity (GPP) from 66 flux sites (Supplementary Fig. 1), (2) 30,369 time-series observations from 4,329 in situ sites since the 1950s, (3) the third generation of the normalized difference vegetation index (NDVI, GIMMS NDVI3g version 1) for 1982–2015 and (4) the NDVI dataset from the MOD13C1 Moderate-Resolution Imaging Spectroradiometer (MODIS) product (collection 6) for 2001–2018.

## Widespread decreases in $P_{freq}$ in northern ecosystems

In the observation records, both winter and spring $P_{freq}$ tended to decrease significantly in the Climatic Research Unit (CRU) gridded time series, the fifth-generation European Centre for Medium-Range Weather Forecasts reanalysis for agriculture and agro-ecological studies (AgERA5) (1982–2018), and the FLUXNET rain gauge data (1989–2014) (Fig. 1a,c). Average $P_{freq}$ and its spatial distribution and temporal pattern were overall consistent for

<footnote n="1">¹Department of Geography, The Ohio State University, Columbus, OH, USA. ²Laboratoire des Sciences du Climat et de l'Environnement, IPSL-LSCE CEA CNRS UVSQ, Gif-sur-Yvette, France. ³CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona, Spain. ⁴CREAF, Cerdanyola del Valles, Barcelona, Spain. ⁵e-mail: wang.12679@osu.edu; liu.738@osu.edu</footnote>
a

<figure><img src="image_3.png" /></figure>

b

<figure><img src="image_4.png" /></figure>

c

<figure><img src="image_5.png" /></figure>

d

<figure><img src="image_6.png" /></figure>

**Fig. 1** | Temporal trends of precipitation frequency (*P*<sub>freq</sub>) in northern ecosystems (>30°N). **a,c**, Trends of winter (December–February; **a**) and spring (March–May; **c**) *P*<sub>freq</sub> anomalies for CRU, AgERA5 (1982–2018, Methods) and FLUXNET (1989–2014) data. **b,d**, Spatial distribution of winter (**b**) and spring (**d**) *P*<sub>freq</sub> trends for average (CRU and AgERA5) data during 1982–2018. P, N and NS indicate the percentages of significantly positive, significantly negative, and non-significant trends, respectively (*P* < 0.05). Grey represents non-significant and none/sparsely vegetated areas.

CRU and AgERA5 (Supplementary Fig. 2), so we used the average (CRU and AgERA5) data as the final *P*<sub>freq</sub>. We found predominantly decreasing trends of winter *P*<sub>freq</sub> (42.7% of the area) and spring *P*<sub>freq</sub> (37.8%) against smaller areas with increasing trends (winter: 9.2%; spring: 7.3%) in northern ecosystems (*P* < 0.05) during 1982–2018 (Fig. 1b,d). Decreasing trends of *P*<sub>freq</sub> were widespread (such as in Siberia and northern Europe) while increasing trends were localized in specific areas such as western Canada and the northern United States.

### Response of LOD to *P*<sub>freq</sub> at different scales

As for trends in LOD, we found that GPP-based LOD of 66 sites significantly advanced and delayed (*P* < 0.05) at nine and two sites, respectively (Supplementary Fig. 3a). Similarly, LOD showed advancing (40.5, 52.2 and 8.6% of the area) and delaying (4.5, 16.1 and 3.5%) trends (*P* < 0.05) for in situ, NDVI3g and MODIS data, respectively (Supplementary Fig. 3b–d). *T*<sub>mean</sub>, *P*<sub>total</sub> and *C*<sub>total</sub> of preseason, the site-dependent period before LOD with the highest absolute partial-correlation coefficient (Methods), have been reported to have larger impacts on LOD than in winter or spring<sup>4,8</sup>. Thus, we applied partial-correlation analyses to investigate the response of LOD to variations of preseason precipitation under three scenarios: (1) LOD versus *P*<sub>total</sub> controlling *T*<sub>mean</sub> and *C*<sub>total</sub> (PARCOR1), (2) LOD versus *P*<sub>total</sub> controlling *T*<sub>mean</sub>, *C*<sub>total</sub> and *P*<sub>freq</sub> (PARCOR2) and (3) LOD versus *P*<sub>freq</sub> controlling *T*<sub>mean</sub>, *C*<sub>total</sub> and *P*<sub>total</sub> (PARCOR3) (Methods and Supplementary Table 1). The partial correlation between anomalies of GPP-based LOD and *P*<sub>total</sub> under PARCOR1 was significantly positive for the 66 sites combined (745 site-year records) (*P* < 0.05), indicative of the strong control of GPP-based LOD variability. Grouping sites into

plant functional types generated similar results, with significant partial correlations for deciduous broadleaf forests (*P* < 0.01) and mixed forests (*P* < 0.05) (Fig. 2a). The overall partial correlation became non-significant, however, after removing the effect of preseason *P*<sub>freq</sub> on GPP-based LOD (PARCOR2) (Fig. 2e). By contrast, positive partial correlations (*P* < 0.001) were overall maintained between anomalies of GPP-based LOD and *P*<sub>freq</sub> under PARCOR3 (Fig. 2i), indicating the importance of *P*<sub>freq</sub> in controlling interannual variability of LOD and the relationship between LOD and *P*<sub>total</sub>.

Analysis of in situ observations of LOD from 4,329 sites for 28 species (total of 30,369 time series) generated similar results. The partial correlation between ground-based LOD and *P*<sub>total</sub> under PARCOR1 was significantly positive (*P* < 0.05) for 14.7% of the time series, nearly twice the number of the significantly negative counterparts (7.3%, Fig. 2b). The total percentages of significant time series decreased to 9.3% under PARCOR2 (Fig. 2f). Yet, 22% of ground-based LOD remained significantly (*P* < 0.05) partially correlated with *P*<sub>freq</sub> under PARCOR3, 64.4% with positive partial correlation (Fig. 2j). Positive-dominant effects of *P*<sub>total</sub> (PARCOR1) on ground-based LOD, especially for typical temperate tree species (*A. hippocastanum* L. and *B. pendula* Roth), agreed with the previous study<sup>14</sup>. Interestingly, we found contrasting effects of *P*<sub>total</sub> (PARCOR1) and *P*<sub>freq</sub> (PARCOR3) on ground-based LOD between temperate tree species (positive-dominant) and meadows (negative-dominant), indicating divergent responses of woody versus herbaceous species to the two precipitation indicators. Sites with significantly negative correlations under PARCOR1 and PARCOR3 were generally located in relatively warm areas (>4 °C) during preseason (Supplementary Fig. 4a,d).
a

<figure><img src="image_7.png" /></figure>

b

<figure><img src="image_8.png" /></figure>

c

<figure><img src="image_9.png" /></figure>

d

<figure><img src="image_10.png" /></figure>

e

<figure><img src="image_11.png" /></figure>

f

<figure><img src="image_12.png" /></figure>

g

<figure><img src="image_13.png" /></figure>

h

<figure><img src="image_14.png" /></figure>

i

<figure><img src="image_15.png" /></figure>

j

<figure><img src="image_16.png" /></figure>

k

<figure><img src="image_17.png" /></figure>

l

<figure><img src="image_18.png" /></figure>

**Fig. 2 | Impact of precipitation on LOD in northern ecosystems (>30° N).** Partial correlations (PARCOR) between LOD and precipitation under three scenarios. **a–d**, PARCOR1. **e–h**, PARCOR2. **i–l**, PARCOR3. **a,e,i**, FLUXNET data. ENF, evergreen needleleaf forests; DBF, deciduous broadleaf forests; MF, mixed forests; GRA, grasslands. **b,f,j**, In situ data. **c,g,k**, NDVI3g data (1982–2015). **d,h,l**, MODIS data (2001–2018). P and N indicate the percentage of significantly positive and negative partial correlations, respectively (*P* < 0.05). Grey represents non-significant and none/sparsely vegetated areas.

Results from the analysis of satellite greenness products agreed with these findings. Partial correlations between NDVI3g-based LOD (1982–2015) and *P*<sub>total</sub> under PARCOR1 were positive (*P* < 0.05) in 22.5% of the area, nearly four times the area with significantly negative correlations (5.8%, Fig. 2c). The total area with significant partial correlation decreased by 49% under PARCOR2 (Fig. 2g). Moreover, 16.7% of the area had significant and positive partial correlations under PARCOR1, more than twice the area with significantly negative correlation for MODIS data (2001–2018) (Fig. 2d). The total areas with significant correlations, however, also decreased by 32% under PARCOR2 (Fig. 2h). As for *P*<sub>freq</sub> effects, 73% and 64% of the area with significant correlation under PARCOR3 were positive for NDVI3g and MODIS data (Fig. 2k,l), respectively. For NDVI3g data, significantly negative correlations under PARCOR1 and PARCOR3 were mainly in warm and dry regions with soil temperatures

>3 °C and soil moisture <0.15 m³ m⁻³ (Supplementary Fig. 4b,e). For MODIS data, negative correlations under PARCOR1 and PARCOR3 were mainly in relatively dry regions (Supplementary Fig. 4c,f). Patterns of PARCOR1 and PARCOR3 were similar in different biomes (Supplementary Fig. 5), and satellite-based LOD for herbaceous biomes (temperate and montane grasslands) and woody biomes showed contrasting responses to *P*<sub>total</sub> and *P*<sub>freq</sub>. To account for the effect of rainfall size in the frequency indicator, we also explored the impact of *P*<sub>freq</sub> for different rainfall event sizes (1 mm d<sup>-1</sup>, 5 mm d<sup>-1</sup> and 10 mm d<sup>-1</sup>) on satellite-based LOD. Two-thirds of the significant correlations between *P*<sub>freq</sub> at 1 mm d<sup>-1</sup> and LOD are positive (*P* < 0.05) under PARCOR3, but this discrepancy became non-existent for *P*<sub>freq</sub> at 5 mm d<sup>-1</sup> and *P*<sub>freq</sub> at 10 mm d<sup>-1</sup> (Supplementary Fig. 6), indicating that the effect of *P*<sub>freq</sub> is controlled by total *P*<sub>freq</sub> rather than by the frequency of large rainfall events. These results suggest that the dominant positive
<figure><img src="image_19.png" /><figcaption>**Fig. 3** | Climatic response to LOD. **a**, Dominant climatic factors for the NDVI3g data (Methods). **b**, Average climate signal, defined as the absolute value of SV. **c–f**, The SVs derived from ridge regression for $T_{\text{mean}}$ (**c**), $P_{\text{freq}}$ (**d**), $P_{\text{total}}$ (**e**) and $C_{\text{total}}$ (**f**). ND and NS indicate no dominant factor and non-significant regression ($P < 0.05$), respectively. Mean indicates the mean value of SV for all significant areas. Positive and negative SV indicate delaying and advancing effects on LOD, respectively. Grey represents non-significant and none/sparsely vegetated areas. The MODIS and in situ results are detailed in Supplementary Figs. 7 and 8, respectively.</figcaption></figure>

partial correlation between LOD and precipitation was mainly influenced by $P_{\text{freq}}$ instead of $P_{\text{total}}$.

## Sensitivity of $P_{\text{freq}}$ to LOD

Analyses of all four independent lines of evidence (carbon-flux measurements, in situ records and data from the NDVI3g and MODIS greenness) confirmed an essential role of $P_{\text{freq}}$ in controlling the effect of precipitation on LOD (previous section). Here we used the climatic signal, calculated as the absolute value of climatic sensitivity (SV, Methods)²⁸, to assess the extent to which climatic factors influence LOD and determine the dominant factor. On the basis of NDVI3g data, we found that, among climatic factors, preseason $P_{\text{freq}}$ dominated over 9.7% of the area, close to $T_{\text{mean}}$ (10.8%), with a larger contribution than $P_{\text{total}}$ and $C_{\text{total}}$ (Fig. 3a,b), suggesting a role of $P_{\text{freq}}$ in explaining LOD variations. Sensitivity analyses indicate that $T_{\text{mean}}$ had a negative-dominant effect on LOD, whereas $P_{\text{freq}}$ had overall positive effects, especially in the high latitudes (Fig. 3c,d). The mean value of sensitivities also indicates the direction and extent to which climatic factors influence LOD. $P_{\text{freq}}$ (0.13) overall had a stronger effect on LOD than $P_{\text{total}}$ (0.02) and $C_{\text{total}}$ (0.02) (Fig. 3d–f). Given the recent widespread decrease in $P_{\text{freq}}$ (Fig. 1), these results also suggest a positive contribution of $P_{\text{freq}}$ to the advance of LOD. Similar results were obtained for MODIS data (Supplementary Fig. 7). For in situ observations, we found similar results that preseason

$P_{\text{freq}}$ showed a stronger influence than $P_{\text{total}}$ and $C_{\text{total}}$ for different species (Supplementary Fig. 8a–f). Interestingly, unlike for temperate tree species, $P_{\text{freq}}$ sensitivity of meadows was negative-dominated (Supplementary Fig. 8g), consistent with the sign of partial correlation between $P_{\text{freq}}$ and LOD (Fig. 2j). Furthermore, LOD in preseasons with lower $P_{\text{freq}}$ exhibits a stronger response to $P_{\text{total}}$ than in preseasons with higher $P_{\text{freq}}$ for in situ and NDVI3g data (Supplementary Fig. 9), indicating a nonlinear response to precipitation controlled by $P_{\text{freq}}$.

## Mechanisms of the effect of $P_{\text{freq}}$

Several mechanisms probably underlay the response of LOD to changes in $P_{\text{freq}}$. First, surface absorbed radiation is directly influenced by $P_{\text{freq}}$, supported by negative-dominant partial correlations between gridded and flux-tower based $P_{\text{freq}}$ and radiation annual variations (Fig. 4a and Supplementary Fig. 10). Nearly 75% of the area with a significant partial correlation between radiation and satellite-based LOD had a negative correlation value (Fig. 4d), indicating that decreases in $P_{\text{freq}}$ associated with less cloudiness, enhance radiation and further lead to earlier LOD. $P_{\text{freq}}$-induced changes in radiation could modulate the heat requirement for leaf unfolding¹⁵, especially when accumulated chilling is not fulfilled. Second, reduced rainfall frequency, accompanied with more clear-sky days and nights, increases the daytime surface solar heating and decreases night-time downward longwave radiation, leading to higher daytime temperature ($T_{\text{max}}$) and lower night-time temperature ($T_{\text{min}}$)²⁹ (Fig. 4b,c). These contrasting effects contribute to earlier LOD with predominantly negative ($T_{\text{max}}$ versus LOD) and positive ($T_{\text{min}}$ versus LOD) partial correlations (Fig. 4e,f), suggesting that widespread decreases in $P_{\text{freq}}$ could concurrently accelerate heat accumulation (at days) and chilling accumulation (at night) before leaf onset. Climatic warming has dual effects on LOD. Specifically, warming could advance LOD, but this effect is counteracted by the reduced chilling during dormancy⁵,⁶. Our results not only support inconsistent responses of LOD to daytime and night-time warming shown in ref. ⁸, but also show a positive contribution of lower $P_{\text{freq}}$ on LOD advancement via synergetic effects on both higher $T_{\text{max}}$ and lower $T_{\text{min}}$.

Notably, almost one-third of significant correlations ($P_{\text{freq}}$ versus LOD) for in situ and satellite data were negative (Fig. 2j-1), meaning that, for example, a decreased $P_{\text{freq}}$ comes with a delayed LOD, requiring alternative explanations than those proposed in the preceding. Grouping correlations into different species (biomes) indicates opposite effects of $P_{\text{freq}}$ on woody (positive-dominant) versus herbaceous (negative-dominant) plants (Fig. 2j and Supplementary Fig. 5c,d). Here we gave a potential mechanism of $P_{\text{freq}}$ effects for grasslands that are located mainly in semiarid regions. Using reanalysis-based soil moisture and a drought indicator (Standardized Precipitation Evapotranspiration Index (SPEI)), we found, after removing the effect of $P_{\text{total}}$, the decreases in $P_{\text{freq}}$ led to lower soil water availability (Supplementary Fig. 11a,c), and concentrated rainfall enhanced water losses from run-off²⁵ (Supplementary Fig. 11b). This drought stress further delayed LOD as shown by predominantly negative correlations (Supplementary Fig. 11d), indicating that decreases in $P_{\text{freq}}$ could aggravate drought stress and delay LOD accordingly in grasslands. This tendency to postpone LOD and associated evapotranspiration could reflect a strategy for herbaceous species³⁰ or some woody species³¹ to adapt to water depletion. Decreased soil moisture might partly reduce nutrient availability (for example, nitrogen) in arid and semiarid regions³²,³³ and further delay LOD¹⁴, requiring additional manipulation experiments. This evidence overall supports our hypothesis that lower $P_{\text{freq}}$ contributes to the advance of LOD in northern ecosystems.

## Modelling and projections of LOD

Most current spring phenological models based solely on daily $T_{\text{mean}}$, such as conventional threshold methods (CT) and growing-degree
<figure><img src="image_20.png" /><figcaption>**Fig. 4** | Mechanisms of the effect of $P_{freq}$ on LOD. **a–f**, Spatial patterns of partial correlations: $P_{freq}$ versus incoming solar radiation ($R_{solar}$) (**a**), $P_{freq}$ versus $T_{max}$ (**b**), $P_{freq}$ versus $T_{min}$ (**c**), $R_{solar}$ versus LOD (**d**), $T_{max}$ versus LOD (**e**) and $T_{min}$ versus LOD (**f**). LOD is derived from NDVI3g data (1982–2015). P and N indicate the percentages of significantly positive and negative partial correlations, respectively ($P < 0.05$). Grey represents non-significant and none/sparsely vegetated areas.</figcaption></figure>

days (GDD), ignore the predictive strength of precipitation in controlling vegetation seasonality⁹. Previous studies have illustrated the importance of precipitation variations in improving the estimation of satellite-based LOD³⁴. Here, we developed a new algorithm called GDDPREC (Methods) for predicting LOD by incorporating information on precipitation ($P_{total}$ and $P_{freq}$) into GDD model, and we compared the performances of CT, GDD and GDDPREC models using both in situ and satellite observations (Fig. 5a–d). The new model (GDDPREC) improved the prediction of frequency of sites/pixels with significant correlation (observational LOD versus predicted LOD, $P < 0.05$), the correlation coefficient ($R$), the root mean square error (RMSE), the corrected Akaike information criterion (AICc, Methods) and the simulation of temporal trends of LOD. A fraction of 82, 61 and 35% of the time series from modelled GDDPREC showed significant positive correlations with observed LOD using in situ, NDVI3g and MODIS data, respectively. These percentages decreased to 37, 39, and 19% for CT models and 66, 51 and 25% for the GDD-only models, respectively (Fig. 5a). Average $R$ indicated 132, 52 and 47% increases versus CT and 32, 23 and 31% increases versus GDD (Fig. 5b). Lower RMSE further confirmed the improvement of LOD modelling by the GDDPREC model (Fig. 5c). The GDDPREC model reduced AICc by 23, 19, and 16% versus CT and 10, 8 and 8% versus GDD using observed LOD from in situ, NDVI3g and MODIS data, respectively (Fig. 5d). In addition, we found a lower absolute difference of LOD regression slope between observed LOD and modelled value from GDDPREC compared with LOD modelled by CT and GDD (Supplementary Fig. 12), indicating the improvement of GDDPREC on predicting the temporal trends of LOD.

Our new model improved the accuracy of LOD prediction, so we applied it to predict future LOD under the representative concentration pathway (RCP) 4.5 and RCP 8.5 future scenarios using temperature

and precipitation bias-corrected model (Supplementary Table 2) projections during 2019–2099 (Fig. 5e–j). Compared with the ensemble mean LOD derived from GDDPREC during 2080–2099, CT advanced LOD estimation in northern Canada and northeastern Asia, with spatially averaged differences of 0.6 and –0.3 d under RCP 4.5 and RCP 8.5, respectively (Fig. 5e,g). Relative to the widely used GDD, the ensemble mean LOD from GDDPREC was predicted to be earlier than currently expected in 62.3% and 68.1% of the area under RCP 4.5 and RCP 8.5 for 2080–2099, respectively (Fig. 5f,h). Grouping the results into biomes yielded overall overestimation of LOD (Fig. 5i). Ensemble mean LOD derived from GDDPREC tended to significantly advance during 2019–2099, with slopes of –0.12 and –0.22 dyr$^{-1}$ under RCP 4.5 and RCP 8.5 ($P < 0.001$), respectively (Fig. 5j). Projections of LOD from individual bias-corrected models showed similar overestimation of LOD (Supplementary Fig. 13), contributing to a negative feedback to climate.

## Conclusion

Our results generally indicate a new but important role of $P_{freq}$ in controlling the effect of precipitation on LOD in northern ecosystems. The synthesis of carbon-flux measurements, in situ records, and data from satellite greenness products suggests that the recent decreases in $P_{freq}$ partially explain the advance of LOD. The significant response of LOD to $P_{total}$, consistent with previous studies¹³,¹⁴, could be considerably negated by controlling the effect of $P_{freq}$, indicating the importance of $P_{freq}$ in the relationship between precipitation and LOD. We further found predominantly positive (nearly two-thirds) partial correlations between $P_{freq}$ and LOD. We considered three mechanisms linking variations in $P_{freq}$ with changes in LOD: (1) lower $P_{freq}$ increases surface absorbed radiation, further advancing LOD; (2) decreases in $P_{freq}$ accompanied with more clear-sky days and nights, result in lower night-time temperature
<figure><img src="image_21.png" /></figure>

**Fig. 5** | Comparison of the three predictive algorithms for modelling and projections of LOD. The three predictive algorithms are the CT, GDD and $GDD_{PREC}$ (Methods). **a–d**, The criteria for evaluating the algorithms include the frequency of sites/areas with significant correlation ($P<0.05$) (**a**), **R** (**b**), the RMSE (**c**) and the AICc (**d**). The legend in **a** applies to all panels. **e–h**, Spatial pattern of LOD differences, $GDD_{PREC}$ – CT (RCP 4.5 (**e**), RCP 8.5 (**g**) $GDD_{PREC}$ – GDD (RCP 4.5 (**f**)), RCP 8.5 (**h**) using bias-corrected multimodel (Supplementary Table 2) projections during 2080–2099. P, N and mean indicate the percentages of positive and negative differences and spatially averaged differences, respectively. **i**, Average differences in LOD (2080–2099) for vegetation types. BF, boreal forests; TD, tundra; TBF, temperate broadleaf forests; TG, temperate grasslands; TCF, temperate coniferous forests; MG, montane grasslands; XS, xeric shrublands; MF, Mediterranean forests (Supplementary Fig. 1). **j**, Temporal trends of predicted LOD (2019–2099) using three algorithms. Shaded areas show the standard deviation of LOD.

and higher daytime temperature. Divergent temperature responses concurrently contribute to the advance of LOD, associated with better fulfilments of both chilling and heat requirements; (3) for herbaceous plants located mainly in semiarid regions, lower $P_{freq}$ could aggravate drought stress and delay LOD accordingly. Our improved model generally projected an earlier LOD than currently expected, advancing nearly twice as fast under RCP 8.5 than under RCP 4.5. The length of future growing seasons and the amount of carbon uptake might be consequently underestimated.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgement, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-022-01285-w.

Received: 22 July 2021; Accepted: 11 January 2022; Published online: 14 February 2022

## References

1. Keenan, T. F. et al. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. *Nat. Clim. Change* **4**, 598–604 (2014).

2. Barichivich, J. et al. Large-scale variations in the vegetation growing season and annual cycle of atmospheric CO₂ at high northern latitudes from 1950 to 2011. *Glob. Change Biol.* **19**, 3167–3183 (2013).

3. Vitasse, Y. et al. Assessing the effects of climate change on the phenology of European temperate trees. *Agr. Forest Meteorol.* **151**, 969–980 (2011).

4. Menzel, A. et al. European phenological response to climate change matches the warming pattern. *Glob. Change Biol.* **12**, 1969–1976 (2006).

5. Fu, Y. H. et al. Declining global warming effects on the phenology of spring leaf unfolding. *Nature* **526**, 104–107 (2015).

6. Wang, H. et al. Overestimation of the effect of climatic warming on spring phenology due to misrepresentation of chilling. *Nat. Commun.* **11**, 4945 (2020).

7. Myneni, R. C. et al. Increased plant growth in the northern high latitudes from 1981 to 1991. *Nature* **386**, 698–702 (1997).

8. Piao, S. et al. Leaf onset in the Northern Hemisphere triggered by daytime temperature. *Nat. Commun.* **6**, 6911 (2015).

9. Richardson, A. D. et al. Influence of spring and autumn phenological transitions on forest ecosystem productivity. *Phil. Trans. R. Soc. B* **365**, 3227–3246 (2010).

10. Piao, S. et al. Plant phenology and global climate change: current progresses and challenges. *Glob. Change Biol.* **25**, 1922–1940 (2019).

11. White, A., Cannell, M. G. R. & Friend, A. D. The high-latitude terrestrial carbon sink: a model analysis. *Glob. Change Biol.* **6**, 227–245 (2000).

12. Piao, S. et al. Net carbon dioxide losses of northern ecosystems in response to autumn warming. *Nature* **451**, 49–53 (2008).

13. Fu, Y. H. et al. Unexpected role of winter precipitation in determining heat requirement for spring vegetation green-up at northern middle and high latitudes. *Glob. Change Biol.* **20**, 3743–3755 (2014).
14. Yun, J. et al. Influence of winter precipitation on spring phenology in boreal forests. *Glob. Change Biol.* **11**, 5176-5187 (2018).

15. Fu, Y. H. et al. Increased heat requirement for leaf flushing in temperate woody species over 1980-2012: effects of chilling, precipitation and insolation. *Glob. Change Biol.* **21**, 2687-2697 (2015).

16. Wipf, S., Stoeckli, V. & Bebi, P. Winter climate change in alpine tundra: plant responses to changes in snow depth and snowmelt timing. *Climatic Change* **94**, 105-121 (2009).

17. Peñuelas, J. et al. Complex spatiotemporal phenological shifts as a response to rainfall changes. *New Phytol.* **161**, 837-846 (2004).

18. Paschalis, A. et al. Rainfall manipulation experiments as simulated by terrestrial biosphere models: where do we stand? *Glob. Change Biol.* **26**, 3336-3355 (2020).

19. Green, J. K. et al. Regionally strong feedbacks between the atmosphere and terrestrial biosphere. *Nat. Geosci.* **10**, 410-414 (2017).

20. Trenberth, K. E., Dai, A., Rasmussen, R. M. & Parsons, D. B. The changing character of precipitation. *Bull. Am. Meteorol. Soc.* **84**, 1205-1217 (2003).

21. Qian, W., Fu, J. & Yan, Z. Decrease of light rain events in summer associated with a warming environment in China during 1961-2005. *Geophys. Res. Lett.* **34**, L11705 (2007).

22. Sun, Y., Solomon, S., Dai, A. & Portmann, R. W. How often will it rain? *J. Clim.* **20**, 4801-4818 (2007).

23. Chou, C. et al. Mechanisms for global warming impacts on precipitation frequency and intensity. *J. Clim.* **13**, 3291-3306 (2012).

24. Held, I. M. & Soden, B. J. Robust responses of the hydrological cycle to global warming. *J. Clim.* **19**, 5686-5699 (2006).

25. Fowler, M. D., Kooperman, G. J., Randerson, J. T. & Pritchard, M. S. The effect of plant physiological responses to rising CO₂ on global streamflow. *Nat. Clim. Change* **9**, 873-879 (2019).

26. Belnap, J., Phillips, S. L. & Miller, M. E. Response of desert biological soil crusts to alterations in precipitation frequency. *Oecologia* **141**, 306-316 (2004).

27. Knapp, A. K. et al. Consequences of more extreme precipitation regimes for terrestrial ecosystems. *Bioscience* **58**, 811-821 (2008).

28. Chen, L. et al. Leaf senescence exhibits stronger climatic responses during warm than during cold autumns. *Nat. Clim. Change* **10**, 777-780 (2020).

29. De Boeck, H. J., Dreesen, F. E., Janssens, I. A. & Nijs, I. Climatic characteristics of heat waves and their simulation in plant experiments. *Glob. Change Biol.* **16**, 1992-2000 (2010).

30. Shen, M. et al. Precipitation impacts on vegetation spring phenology on the Tibetan Plateau. *Glob. Change Biol.* **21**, 3647-3656 (2015).

31. Peaucelle, M. et al. Spatial variance of spring phenology in temperate deciduous forests is constrained by background climatic conditions. *Nat. Commun.* **10**, 5388 (2019).

32. Estiarte, M. & Peñuelas, J. Alteration of the phenology of leaf senescence and fall in winter deciduous species by climate change: effects on nutrient proficiency. *Glob. Change Biol.* **21**, 1005-1017 (2015).

33. Austin, A. T. et al. Water pulses and biogeochemical cycles in arid and semiarid ecosystems. *Oecologia* **141**, 221-235 (2004).

34. White, M. A., Thornton, P. E. & Running, S. W. A continental phenology model for monitoring vegetation responses to interannual climatic variability. *Glob. Biogeochem. Cycles* **11**, 217-234 (1997).

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2022
Methods

**In situ observations.** We applied three independent in situ datasets for ground-based LOD (leaf unfolding date, LUD) (>30° N):

(1) The Pan European Phenology Project³⁵ (PEP725, http://www.pep725.eu/), which provides an open-access and long-term (since 1868) phenological database for 19,608 sites and 78 species across 25 European countries

(2) The Chinese Phenological Observation Network³⁶ (CPON), which has compiled phenological observations since 1963 for 112 species and 145 sites across China

(3) The USA National Phenology Network³⁷ (NPN, https://www.usanpn.org/), which has received contributions from many citizen scientists using a standardized protocol for observing plant phenology across the United States

The definition of spring LUD differs among the three datasets. PEP725, CPON and NPN define LUD as the date of the first visible foliar stalk for tree species (BBCH code 11) and 25% green in spring for meadow (BBCH code 101), 50% full foliar expansion and the timing of the first bud break, respectively. To identify and remove potential outliers, we applied the median absolute deviation (MAD) method, which is more resilient to outliers in a dataset than the standard deviation. In our case, MAD of LUD dataset (LUD₁, LUD₂, ..., LUDₜ) can be expressed as:

$$MAD = \text{median}(|\text{LUD}_t - \text{median} (\text{LUD})|) \quad (1)$$

For each site, any data record with more than 2.5 times MAD is considered an outlier. We also excluded all LUD records that were shorter than 15 years. In this way, we used a total of 30,369 time series from 4,329 sites and 28 species for 1951–2018. The distribution and descriptions of the in situ sites are detailed in Supplementary Fig. 1 and Table 3.

**Carbon-flux phenology.** We used eddy-covariance flux measurements to determine the GPP-based LOD (the start of growing season, SOS). After removing sites with insufficient observations (<5 yr), we applied all 66 available flux sites (Supplementary Fig. 1 and Table 4) with a total of 745 year-site records of daily GPP from the FLUXNET database (https://fluxnet.org/). We applied a site-based relative threshold of 10% of the annual maximum GPP to determine SOS³⁰. The choice of relative threshold does not affect the interannual variability of SOS, but higher or lower thresholds will lead to later or earlier mean SOS, respectively¹. We thus utilized yearly anomalies of SOS from all sites for the same plant function type to analyse the responses of SOS to precipitation at the plant-type level.

**Satellite-based phenology.** Two independent satellite greenness products were applied to determine the satellite-based LOD (vegetation green-up date, VGD). GIMMS NDVI3g v1 data (1982–2015) were derived from the measurements of advanced very high resolution radiometer having a spatial resolution of 1/12° and a temporal resolution of 15 days. Terra MODIS NDVI data (2001–2018) were derived from the 16-day MOD13C1 composite product³⁹ (collection 6) with a spatial resolution of 0.05°.

To exclude snow effects, we substituted all contaminated NDVI by the mean of snow-free NDVI values in winter (December–February) of all years⁴⁰. A modified Savitzky–Golay filter was then applied to remove the abnormal values and reconstruct NDVI time series⁴¹. In addition, we eliminated areas with sparse vegetation by removing areas with a mean annual NDVI <0.1 (ref. ³²). We applied two methods to calculate VGD to minimize the uncertainty from a single method, the dynamic-threshold approach and the double-logistic function⁴³.

We calculated NDVI ratios annually for each pixel as:

$$\text{NDVI}_{\text{ratio}} = \frac{\text{NDVI} - \text{NDVI}_{\text{min}}}{\text{NDVI}_{\text{max}} - \text{NDVI}_{\text{min}}} \quad (2)$$

where NDVI, NDVI<sub>min</sub>, and NDVI<sub>max</sub> are the daily NDVI and the annual minimum and maximum of the NDVI curve, respectively. Spring VGD was defined as the day of the year when the NDVI<sub>ratio</sub> increased to 0.5 (ref. ³⁴).

We divided the annual NDVI curve into two sections using the maximum NDVI and applied a piecewise logistic function to fit each section for each area⁴⁴.

$$y(t) = a_1 + (a_2 - a_7)t \left[ \frac{1}{1 + e^{(a_3 - t)/a_4}} - \frac{1}{1 + e^{(a_5 - t)/a_6}} \right] \quad (3)$$

where *t* is time in days, *y*(*t*) is the NDVI at time *t* and *a*<sub>1</sub>–*a*<sub>7</sub> are fitting parameters: *a*<sub>1</sub> is the background NDVI; *a*<sub>2</sub> is the difference between the background and the amplitude of the late summer and autumn plateau, both in NDVI units; *a*<sub>3</sub> and *a*<sub>5</sub> are the midpoints in the days of the year of the transitions for green-up and senescence/abscission, respectively; *a*<sub>4</sub> and *a*<sub>6</sub> are the transition curvature parameters (normalized slope coefficients); and *a*<sub>7</sub> is the summer green-down parameter. Spring VGD was defined as the time when the rate of change in curvature reached its first local maximum in spring.

These two methods produce similar results³⁴, so we determined average VGD from the dynamic-threshold approach and double-logistic function as the final satellite-based LOD. To exclude the impact of human activity on agricultural ecosystems, we removed all cropland areas using the MCD12Q1

MODIS land-cover product (collection 6). We then utilized the borders of the biomes⁴⁵ to conduct the analyses for different vegetation types (Supplementary Fig. 1). Some caution is needed when interpreting the results for heterogeneous pixels within different biomes. It also should be noted that there could be some biases between ground-, GPP- and satellite-based LOD, especially regarding the photosynthesis processes and greenness changes. To minimize this effect, we conducted independent analyses for different datasets (carbon-flux measurements, in situ records and data from two satellite greenness products) instead of directly integrating or comparing these datasets.

**Climatic data.** We derived two independent datasets of precipitation frequency (*P*<sub>freq</sub>, number of rainy days per month) from (1) the CRU time series⁴⁶ (CRU-TS 4.03) at a spatial resolution of 0.5° (https://rudua.uea.ac.uk/cru/data/hrg/), which is interpolated by massive climatic stations, and (2) the AgERA5 at a spatial resolution of 0.1° (https://cds.climate.copernicus.eu). CRU provides a monthly climatological variable of the number of rainy days, defined as the number of rainy days with ≥0.1 mm of precipitation²²,³³,⁴⁷. We extracted AgERA5-based monthly numbers of rainy days using daily AgERA5 precipitation (≥0.1 mm). We noticed that multiyear averages and trends of *P*<sub>freq</sub> from CRU and AgERA5 were very similar (Supplementary Fig. 2), so we calculated the average *P*<sub>freq</sub> and *P*<sub>total</sub> (mm month⁻¹) datasets for CRU and AgERA5 as final *P*<sub>freq</sub> and *P*<sub>total</sub> for 1982–2018 to reduce the uncertainty from a single dataset. Monthly *P*<sub>freq</sub> and *P*<sub>total</sub> during 1950–1982, monthly surface *T*<sub>mean</sub> (°C) and *C*<sub>total</sub> (‰, a proxy of solar radiation) for 1951–2018 and monthly *T*<sub>max</sub> (°C) and *T*<sub>min</sub> (°C) for 1982–2015 at a spatial resolution of 0.5° were obtained from CRU. For the flux sites, we directly utilized monthly *T*<sub>mean</sub>, incoming shortwave radiation (W m⁻²), *P*<sub>total</sub> and *P*<sub>freq</sub> (number of rainy days with ≥0.1 mm of precipitation) measured by flux towers. For the LOD models, we used daily *T*<sub>mean</sub> (the average of *T*<sub>max</sub> and *T*<sub>min</sub>) and *P*<sub>total</sub> at spatial resolutions of 0.5° from the Climate Prediction Center, provided by the National Oceanic and Atmospheric Administration/Oceanic and Atmospheric Research/Earth System Research Laboratories Physical Sciences Laboratory (https://psl.noaa.gov/). For projections of future LOD under two climatic scenarios (RCP 4.5 and RCP 8.5), we used daily *T*<sub>mean</sub> and *P*<sub>total</sub> (with a spatial resolution of 0.5°×0.5°) simulated by four bias-corrected models from the Inter-Sectoral Impact Model Intercomparison Project ¹⁸ (Supplementary Table 2).

Monthly run-off data for 1982–2015 were derived from TerraClimate⁴⁹, a dataset of monthly climate for global terrestrial surfaces at a spatial resolution of 1/24°. We utilized the monthly SPEI (3-month scalar) for 1982–2015 at a spatial resolution of 0.5°, calculated by the difference between precipitation and potential evapotranspiration from the SPEI base v.2.5 at Consejo Superior de Investigaciones Científicas⁵⁰. Volumetric soil water (a proxy for soil moisture, m³ m⁻³) was derived from ERA5-Land monthly average data. We calculated the average volumetric soil water of the top two layers (0–7 cm, 9–28 cm) as the final monthly soil moisture for mechanistic analyses of herbaceous plants.

**Analyses.** We applied the Theil–Sen slope estimator, a non-parametric and median-based slope estimator, to analyse the past and projected temporal trends of LOD for the ground and satellite observations. The trends were evaluated using the Mann–Kendall trend test at a significance level of 0.05.

*T*<sub>mean</sub>, *P*<sub>total</sub> and *C*<sub>total</sub> jointly control LOD so that a simple linear-correlation analysis would have uncertainties of factor-combined effect. For example, *T*<sub>mean</sub> is numerically related to both LOD and *P*<sub>total</sub>, violating the independence of variables in correlation analyses. We thus applied partial-correlation analysis to explore and explain the impact of *P*<sub>freq</sub> on LOD. The partial-correlation analysis was categorized into three scenarios: (1) partial correlation between LOD and *P*<sub>total</sub>, removing the effects of *T*<sub>mean</sub> and *C*<sub>total</sub> (PARCOR1); (2) partial correlation between LOD and *P*<sub>total</sub>, removing the effects of *T*<sub>mean</sub>, *C*<sub>total</sub> and *P*<sub>freq</sub> (PARCOR2); and (3) partial correlation between LOD and *P*<sub>freq</sub>, removing the effects of *T*<sub>mean</sub>, *C*<sub>total</sub> and *P*<sub>total</sub> (PARCOR3) (Supplementary Table 1). Significance was set at *P* < 0.05, with an *R* threshold of ±0.355 for a 34 yr analysis (NDVI3g, 1982–2015) and ±0.514 for an 18 yr analysis (MODIS, 2001–2018). Preseason forcings predicted LOD better than winter or spring climatic forcing alone; the optimal preseason length differs among species and locations. The preseason period was defined as the period with one-month steps until December of the previous year before the month of multiyear mean LOD. During preseason, the absolute partial-correlation coefficient between LOD and climatic factor (for example, *P*<sub>freq</sub>) should be the highest compared with other periods⁵².

To avoid potential multicollinearity between climatic factors, we applied ridge regression that adds a penalty parameter to reduce the variance of the regression coefficient to determine climatic sensitivities. The response variable was LOD, and the predictors were preseason climatic factors. We used normalized anomalies of climatic factors and LOD as regression inputs, and regression coefficients were determined as climatic sensitivities (SVs), including SV-*T*<sub>mean</sub>, SV-*P*<sub>freq</sub>, SV-*P*<sub>total</sub> and SV-*C*<sub>total</sub>. To directly compare the effect of different climatic factors on LOD, we calculated the absolute value of regression coefficients as climatic signals²⁸, indicating the extent to which climatic factors influence leaf unfolding without considering the direction of the effect (delay, advance). For each pixel, we defined the dominant factor as the factor with the highest climatic signal that is greater than the sum of climatic signals of the other three factors.
To evaluate the LOD models, we calculated the frequency of sites/pixels with significant correlations, R, RMSE, AICc and temporal trends of LOD for CT, GDD and GDD<sub>PREC</sub>, respectively. In our case, the sample size (time series for a site or pixel) was small, so we used AICc to address the potential overfitting of AIC. AICc of the model is:

$$AIC = \frac{2k - 2\hat{L}}{n} \qquad (4)$$

$$\text{where } \hat{L} = -\frac{n}{2}(1 + \ln(2\pi) + \ln\left(\frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{n}\right)), \qquad (5)$$

$$\text{so AICc} = \text{AIC} + \frac{2k^2 + 2k}{n - k - 1} \qquad (6)$$

where k is the number of parameters in the model, n is the sample size, L̂ is the log of the maximized value of the likelihood function for the model, yᵢ is the LOD predicted by the model for year i and ŷᵢ is the estimated LOD based on yᵢ.

**Models for predicting LOD.** Most phenological modules in current ecosystem models are based solely on T<sub>mean</sub>. Previous studies have applied temperature-threshold models (for example, T<sub>mean</sub> > 5°C for five consecutive days<sup>51,52)</sup> to estimate plant spring phenology. GDD models are widely used to estimate past and future spring phenology<sup>53</sup>. Considering the potential impacts of precipitation on LOD, we incorporated precipitation (P<sub>total</sub> and P<sub>freq</sub>) into one of the GDD models (GDD<sub>PREC</sub>) and compared GDD<sub>PREC</sub> with the currently applied CT and GDD model.

We compared the three algorithms (CT, GDD and GDD<sub>PREC</sub>) for LOD estimation using in situ and satellite observations. We calculated the average daily T<sub>mean</sub> of five consecutive days before LOD each year. We then set the multiyear mean as the threshold temperature (T<sub>THOLD</sub>) to predict CT-based LOD. If T<sub>mean</sub> was higher than T<sub>THOLD</sub> for five consecutive days from 1 December of the previous year, the first date was determined as CT-based LOD.

The GDD model was calculated as:

$$\text{GDD}(d) = \max(T_{\text{mean}}(d) - T_b, 0) \qquad (7)$$

$$\text{GDD}_{\text{threshold}} = \sum_{d=d_0}^{\text{LOD}} \text{GDD}(d) \qquad (8)$$

where GDD(d) is the growing degree on date d, T<sub>b</sub> is the base temperature, set as 0°C (5 and 10°C provided similar results in this study), T<sub>mean</sub>(d) is the daily mean temperature on date d, GDD<sub>threshold</sub> is the accumulated growing degree from d<sub>0</sub> to LOD required for leaf unfolding and d<sub>0</sub> is the first day of accumulation, set as 1 December of the previous year. GDD-based LOD was defined as the date that GDD(d) first exceeded the multiyear mean GDD<sub>threshold</sub>.

We incorporated P<sub>total</sub> and P<sub>freq</sub> into the GDD model to predict LOD. We first calculated the multiyear average intensity of precipitation as:

$$\text{AIP} = \text{mean} \left( \frac{\sum_{d=d_0}^{\text{LOD}} P_{\text{total}}(d)}{\sum_{d=d_0}^{\text{LOD}} P_{\text{freq}}(d)} \right) \qquad (9)$$

$$\text{GDD}_{\text{pr}}(d) = \max \left( T_{\text{mean}}(d) + k \times \frac{P_{\text{total}}(d)}{\text{AIP}} - T_b, 0 \right) \qquad (10)$$

where AIP represents the multiyear average intensity of precipitation (mm d⁻¹), d<sub>0</sub> is set as 1 December of the previous year and k is a weighted factor ranging from -15 to 15 with steps of 0.1. The effect of precipitation on LOD prediction is jointly controlled by k, P<sub>total</sub> and P<sub>freq</sub>. Intensive precipitation strongly affected GDD<sub>PREC</sub> (P<sub>total</sub>(d) > 1). If P<sub>total</sub> on date d was 0 or k was 0, the accumulated growing degree was solely dependent on T<sub>mean</sub>.

We selected the optimal parameters for GDD<sub>PREC</sub> by comparing the RMSEs between the modelled and observed LOD. k with the lowest RMSE was determined as the final weighted factor. We used the map of k and GDD<sub>threshold</sub> based on GDD<sub>PREC</sub> for 1982–2015 as empirical input data to predict LOD for 2019–2099 (Supplementary Fig. 14).

**Reporting Summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

The in situ phenological data can be accessed from http://www.pep725.eu/ and https://www.usanpn.org/. The flux datasets can be accessed from https://fluxnet.org/. The MODIS NDVI datasets can be accessed from https://modis.gsfc.nasa.gov/data/dataprod/mod13.php. The CRU TS4.00 datasets can be accessed from https://crudata.uea.ac.uk/cru/data/hrg/. The AgERA5 data can be accessed from https://

cds.climate.copernicus.eu. The TerraClimate data can be accessed from http://www.climatologylab.org/terraclimate.html. The CPC datasets can be accessed from https://psl.noaa.gov/. The data for future climates (2019–2099) are available at https://esg.pik-potsdam.de/search/isimip/.

## Code availability

The codes used for data analysis in this study are available on Zenodo at https://doi.org/10.5281/zenodo.5801049.

## References

35. Templ, B. et al. Pan European Phenological database (PEP725): a single point of access for European data. *Int. J. Biometeorol.* **62**, 1109–1113 (2018).

36. Ge, Q., Wang, H., Rutishauser, T. & Dai, J. Phenological response to climate change in China: a meta-analysis. *Glob. Change Biol.* **21**, 265–274 (2015).

37. Schwartz, M. D., Betancourt, J. L. & Weltzin, J. F. From Caprio's lilacs to the USA National Phenology Network. *Front. Ecol. Environ.* **10**, 324–327 (2012).

38. Wu, C. et al. Interannual variability of net ecosystem productivity in forests is explained by carbon flux phenology in autumn. *Glob. Ecol. Biogeogr.* **22**, 994–1006 (2013).

39. Zhang, X. et al. Monitoring vegetation phenology using MODIS. *Remote Sens. Environ.* **84**, 471–475 (2003).

40. Shen, M. et al. Can changes in autumn phenology facilitate earlier green-up date of northern vegetation? *Agr. Forest Meteorol.* **291**, 108077 (2020).

41. Chen, J. et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. *Remote Sens. Environ.* **91**, 332–344 (2004).

42. Shen, M. et al. Increasing altitudinal gradient of spring vegetation phenology during the last decade on the Qinghai–Tibetan Plateau. *Agr. Forest Meteorol.* **189**, 71–80 (2014).

43. Wu, C. et al. Widespread decline in winds delayed autumn foliar senescence over high latitudes. *Proc. Natl Acad. Sci. USA* **118**, e2015821118 (2021).

44. Elmore, A. J. et al. Landscape controls on the timing of spring, autumn, and growing season length in mid-Atlantic forests. *Glob. Change Biol.* **18**, 656–674 (2012).

45. Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on Earth. *Bioscience* **51**, 933–938 (2001).

46. Harris, I., Jones, P. D., Osborn, T. J. & Lister, D. H. Updated high-resolution grids of monthly climatic observations—the CRU TS3.10 dataset. *Int. J. Climatol.* **34**, 623–642 (2014).

47. New, M., Hulme, M. & Jones, P. D. Representing twentieth-century space–time climate variability. Part I: development of a 1961–90 mean monthly terrestrial climatology. *J. Clim.* **12**, 829–856 (1999).

48. Hempel, S., Frieler, K., Warszawski, L., Schewe, J. & Piontek, F. A trend-preserving bias correction—the ISI-MIP approach. *Earth Syst. Dynam.* **4**, 219–236 (2013).

49. Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. *Sci. Data* **5**, 170191 (2018).

50. Vicenteserrano, S. M., Begueria, S. & López-Moreno, J. I. A multiscalar drought index sensitive to global warming: the Standardized Precipitation Evapotranspiration Index. *J. Clim.* **23**, 1696–1718 (2010).

51. Barr, A. G. et al. Inter-annual variability in the leaf area index of a boreal aspen–hazelnut forest in relation to net ecosystem production. *Agr. Forest Meteorol.* **126**, 237–255 (2004).

52. Chen, J., Chen, W., Liu, J., Cihlar, J. & Gray, S. Annual carbon balance of Canada's forests during 1895–1996. *Glob. Biogeochem. Cycles* **14**, 839–849 (2000).

53. Krinner, G. et al. A dynamic global vegetation model for studies of the coupled atmosphere–biosphere system. *Glob. Biogeochem. Cycles* **19**, GB1015 (2005).

## Acknowledgements

This study was funded by the National Science Foundation (no. 1724786). We appreciate principal investigators of flux sites for providing their valuable data for our analyses. We acknowledge all members of the PEP725, CPON and NPN networks for collecting and providing the phenological data. P.C. acknowledges support from the French state aid managed by the ANR under the Investissements d'avenir programme with the reference ANR-16-CONV-0003. J.P. was funded by the Spanish government grant PID2019-110521GB-100, the Fundación Ramón Areces grant ELEMENTAL-CLIMATE and the Catalan government grant SGR2017-1005.

## Author contributions

J.W. and D.L. designed the research. J.W. performed research and analysed data. J.W. wrote the first draft of the manuscript. D.L., P.C. and J.P. substantially revised the manuscript with intensive suggestions.
**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41558-022-01285-w.

**Correspondence and requests for materials** should be addressed to Jian Wang or Desheng Liu.

**Peer review information** *Nature Climate Change* thanks Yongguang Zhang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.
# Reporting Summary

Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see our [Editorial Policies](#) and the [Editorial Policy Checklist](#).

## Statistics

For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a | **Confirmed**

☐ ☑ The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement

☐ ☑ A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

☑ ☐ The statistical test(s) used AND whether they are one- or two-sided

*Only common tests should be described solely by name; describe more complex techniques in the Methods section.*

☑ ☐ A description of all covariates tested

☑ ☐ A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons

☐ ☑ A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

☐ ☑ For null hypothesis testing, the test statistic (e.g. *F*, *t*, *r*) with confidence intervals, effect sizes, degrees of freedom and *P* value noted
*Give *P* values as exact values whenever suitable.*

☑ ☐ For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings

☑ ☐ For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes

☑ ☐ Estimates of effect sizes (e.g. Cohen's *d*, Pearson's *r*), indicating how they were calculated

*Our web collection on statistics for biologists contains articles on many of the points above.*

## Software and code

Policy information about [availability of computer code](#)

Data collection | No software was used.

Data analysis | Arcgis 10.3 was used to produce figures and raster data processing. IDL from ENVI 5.1 was used to calculate satellite-based phenology, partial correlation analysis, and modification of GDD model. Matlab R2016a was used to apply ridge regression for determining climatic sensitivities.

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research [guidelines for submitting code & software](#) for further information.

## Data

Policy information about [availability of data](#)

All manuscripts must include a [data availability statement](#). This statement should provide the following information, where applicable:

- Accession codes, unique identifiers, or web links for publicly available datasets

- A list of figures that have associated raw data

- A description of any restrictions on data availability

The in situ phenological data can be accessed from http://www.pep725.eu/ and https://www.usanpn.org/. The flux data sets can be accessed from https://fluxnet.org/. The data from GIMMS NDVI3g version1 can be accessed from https://ecocast.arc.nasa.gov/data/pub/gimms/. The MODIS NDVI data sets can be accessed from https://modis.gsfc.nasa.gov/data/dataprod/mod13.php. The CRU TS4.00 data sets can be accessed from https://sites.uea.ac.uk/. The AgERA5 data can be accessed from https://cds.climate.copernicus.eu. The TerraClimate data can be accessed from http://www.climatologylab.org/terraclimate.html. The CPC data sets can be accessed from https://psl.noaa.gov/. The data for future climates (2019-2099) is available at https://esg.pik-potsdam.de/search/isimip/.
# Field-specific reporting

Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

<figure><img src="image_22.png" /></figure>

Life sciences

<figure><img src="image_23.png" /></figure>

Behavioural & social sciences

<figure><img src="image_24.png" /></figure>

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see [nature.com/documents/nr-reporting-summary-flat.pdf](https://nature.com/documents/nr-reporting-summary-flat.pdf)

# Ecological, evolutionary & environmental sciences study design

All studies must disclose on these points even when the disclosure is negative.

**Study description**
This research studied the effects of precipitation frequency on spring leaf unfolding date and underlying mechanisms. As a novel aspect, we analyzed and compared the response of spring leaf unfolding to precipitation amount, precipitation frequency, temperature, and solar radiation using eddy covariance flux data, ground observations, and two satellite greenness products. We also improved current GDD model by incorporating precipitation, and predicted future spring leaf unfolding date under RCP4.5 and RCP8.5.

**Research sample**
1. 30,369 time series of phenological observations at 4329 sites since 1950s
2. 745 site-year records of flux measurements at 66 sites
3. GIMMS NDVI 3g data (1982-2015)
4. MODIS NDVI data (2001-2018)

**Sampling strategy**
Because the analyzed data in our study were obtained from open-access database instead of designed experiments, sampling is not applicable to our study.

**Data collection**
Three independent in situ data sets for spring leaf unfolding date were collected from: 1) The Pan European Phenology Project (PEP725, http://www.pep725.eu/). 2) The Chinese Phenological Observation Network. 3) The USA National Phenology Network (NPN, https://www.usanpn.org/).
Daily gross primary productivity (GPP) from the FLUXNET database (www.fluxnet.fluxdata.org).
Two satellite greenness products were collected from GIMMS NDVI 3g and MODIS (MOD13C1).
Climate data were collected from Climatic Research Unit Time Series (CRU-TS 4.03) and the fifth generation ECMWF re-analysis for agriculture and agro-ecological studies (AgERA5)

**Timing and spatial scale**
Between 1951-2018.
Mid-high northern ecosystem (>30 N)

**Data exclusions**
Median absolute deviation (MAD) was used to exclude outliers from the in situ datasets.

**Reproducibility**
Both the phenology, climate data and models are open access. The findings of our study can be reproduced using the statistical methods shown in the manuscript.

**Randomization**
Because the analyzed data in our study were obtained from open-access database instead of designed experiments, randomization is not applicable to our study.

**Blinding**
Because the analyzed data in our study were obtained from open-access database instead of designed experiments, blinding is not applicable to our study.

Did the study involve field work?

<figure><img src="image_25.png" /></figure>

Yes

<figure><img src="image_26.png" /></figure>

No

# Reporting for specific materials, systems and methods

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

## Materials & experimental systems

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☑</td><td>Antibodies</td></tr><tr><td>☑</td><td>Eukaryotic cell lines</td></tr><tr><td>☑</td><td>Palaeontology and archaeology</td></tr><tr><td>☑</td><td>Animals and other organisms</td></tr><tr><td>☑</td><td>Human research participants</td></tr><tr><td>☑</td><td>Clinical data</td></tr><tr><td>☑</td><td>Dual use research of concern</td></tr></tbody></table>

## Methods

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☑</td><td>ChIP-seq</td></tr><tr><td>☑</td><td>Flow cytometry</td></tr><tr><td>☑</td><td>MRI-based neuroimaging</td></tr></tbody></table>