Alignment of tree phenology and climate seasonality influences the runoff response to forest cover loss

Trees shape the critical zone and modulate terrestrial water storage yet observed streamflow responses to forest cover change vary. Differences in catchment area, soil water storage, management practices, tree species, and climate are among the many explanations proposed for heterogeneous hydrologic responses. We addressed evidence for the hypothesis that mean annual temperature (MAT) and the phase shift between precipitation and enhanced vegetation index (EVI) peaks, θ, explain a significant amount of the variation in hydrologic response to forest cover loss. We selected 50 catchments with daily streamflow records spanning eight nations and seven climate regions. Categorical clustering of catchments was performed with MAT, θ, minimum EVI, catchment area, and percentage forest loss. Similar storm event runoff ratio responses to deforestation were best clustered by MAT and θ. High MAT tropical monsoonal catchments (Brazil, Myanmar, and Liberia) exhibited minimal evidence of increasing runoff ratios (increases observed in 9% of catchments). Low MAT subarctic, cold semi-arid, and humid continental catchments (US, Canada, and Estonia) showed consistent runoff increases around the time of snowmelt (94%). The deforestation runoff responses of temperate and subtropical catchments with Mediterranean, humid, and oceanic climates depended strongly on θ. We observe increased runoff following forest loss in a majority of catchments (90%) where precipitation peaks followed peak growing season (max EVI) (US). In contrast, where precipitation peaks preceded the growing season (South Africa and Australia) there was less evidence of increased runoff (25% of catchments). This research supports the strategic implementation of native forest conservation or restoration for simultaneously mitigating the effects of global climate change and regional or local surface runoff.


Introduction
Many nations face complex choices between forest removal to support local economies and food security versus forest preservation to provide a host of ecosystem services including flood control (Ellison et al 2017), water quality (Miura et al 2015), human health (e.g. Mapulanga and Naito 2019), and reduction of global atmospheric carbon (Griscom et al 2017). Forest preservation competes with economic growth in developing (e.g. Ingalls and Dwyer 2016, Vijay et al 2018, Myers et al 2018 and developed nations (Jiang et al 2016, Collentine and Futter 2018), thus the benefits of forest cover must be carefully defined. Observing and explaining catchment-scale hydrologic consequences of forest cover change has remained a challenge (Rogger et al 2017, Mcdonnell et al 2018 and has been highlighted as one of the most critical knowledge gaps facing hydrologic research (Blöschl et al 2019).
Prior studies of forest cover change fall broadly into two categories-those that tracked changes in runoff (i.e. the partitioning of rainfall between immediate streamflow or catchment-stored water) (e.g. Hewlett and Bosch 1984, Trimble et al 1987, Brown et al 2013  . Recent global meta-analyses of streamflow responses to forest cover change indicate that forest loss suppresses evapotranspiration and increases streamflow, though there is considerable noise in these signals indicating heterogeneous hydrologic responses among catchments , Wei et al 2018. Regional studies demonstrate the potential for local hydrologic variations in streamflow-transpiration partitioning (e.g. Peña-Arancibia et al 2019, Goeking and Tarbaton 2020). Together these compilations of forest change research illustrate global heterogeneity in hydrologic responses to forest cover change that has not been adequately explained.
Catchment studies focusing on runoff generation suggest that much of the variability in runoff responses, and possibly water yield, may be due to differences in the timing of annual precipitation extremes, catchment wetness, and the phenology of local vegetation. Ye et al (2017) demonstrated that regions of the contiguous US (CONUS) with energy and water cycle phase shifts show a strong influence of antecedent catchment saturation on the annual maximum flood, whereas catchments with highly synchronized cycles were less sensitive to saturation. Hydrologic studies of forest cover loss across the tropics ( (Blöschl et al 2015, Carrick et al 2019 European catchments have indicated that regional forest cover change could have minimal impact on runoff dynamics due to low water use by vegetation during cool seasons.
Taken together these studies led us to hypothesize that the runoff response to forest cover change may be related to the seasonality of catchment saturation, driven by precipitation and tree phenology. Given the importance of developing a mechanistic understanding of catchment runoff responses to forest cover change for forest management planning (Rogger et al 2017, Ellison et al 2017, we investigate the following questions in this research: • Are seasonal surface runoff responses to forest cover change related to temporal alignments of precipitation peaks and tree phenology?

Hydrologic responses to deforestation
Baseflow-runoff separation was performed for daily stream discharge records with the algorithm of Nathan and Mcmahon (1990)  , where Q SURF is surface runoff and GP is gross precipitation, for each event exceeding 5 mm of precipitation (discussed in section S1, figure S1), generating at least 1 mm of surface runoff. A sensitivity analysis of event rainfall thresholds (5, 10, and 20 mm) indicated that results were not sensitive to threshold choices above 5 mm. Rainfall events were separated with a 48-h minimum inter-event time. We divided each streamflow record into two periods, years preceding and following the largest forest cover loss event between 2000 and 2018 as defined in Hanson et al (2013) (table S1). We tested for changes in paired (preand post-deforestation) storm event runoff ratios by month with a one-sided Mann-Whitney U-test. We defined results as significant when exceeding the α of 0.1 threshold. We grouped catchments into two sets of responses: catchments that experienced (1) a significant storm event runoff ratio increase in any month following forest loss, and (2) no significant increases.
The influence of forest loss on average annual water yield was analyzed in the context of the Budyko model (Budyko 1958), a prediction of catchment partitioning of gross precipitation (GP) between discharge and Actual EvapoTranspiration (AET) (i.e. the Evaporative Index ) as a function of PET ) (equation (1)). (1) We compare pre-and post-deforestation water partitioning with the estimated Budkyo parameter, ω, fit to each cluster (as defined by mean annual temperature (MAT) and θ) with standard least squares.

Hierarchical clustering of catchments by climate, tree phenology, area, and forest loss
Given the size of our dataset, we used hierarchical clustering to group (Akman et al 2019) significant and non-significant runoff responses. For each catchment, we characterized local tree phenology in years before deforestation using the median month of peak EVI. We estimated the timing of annual precipitation extremes by computing the frequency of overthreshold daily rainfall events (m = 5 mm) within bi-weekly intervals in each year. Sensitivity analysis indicated results were not sensitive to the choice window size (1, 2, or 4 weeks). The timing of intense precipitation was estimated as the 2-week window with the most over-threshold events across all years. Phase shift was calculated between precipitation extremes (defined by the starting Julian date of the 2-week period) and EVI (defined as the Julian date of the first day of the month), θ, for each catchment as the Julian date of peak precipitation minus the Julian date of peak EVI.
We compared the predictive power of θ, MAT, median annual minimum EVI, catchment area, and percentage forest loss to explain increases in RR following forest loss. We constructed seven hierarchical clustering approaches, which are described in detail in the supplemental material. Distance between catchments in each centered and scaled multivariate space was defined by the Euclidean norm. The optimal number of clusters was determined with the Calinski-Harabasz Index (Calinski and Harabasz 1974). The ability of cluster models to identify a significant increase in RR was assessed with both the Purity Index (PI) and Normalized Mutual Information (NMI). PI and NMI are objective measures of how effectively clustering can group variables by hidden classes (i.e. significant increases in RR vs. no RR response) (Manning et al 2008). Both PI and NMI range from 0 to 1, where 0 indicates poor class characterization by clusters and 1 indicates perfect characterization. PI and NMI do not require a priori definition of which classes belong in groups. PI compares values within clusters to the mode value. NMI is similar to PI but penalizes clustering as the number of clusters increases.

Runoff responses to deforestation
Shifts in runoff ratios following deforestation varied across countries and climate regions (figures 1 and S1). High latitude catchments in Estonia (figure 1(a)) Canada, and the US (figure 1(b)) exhibited consistent significant runoff ratio increases (increases observed in 16 of 17 of catchments, 94%). In contrast, tropical catchments in Brazil (figure 1(e)), Liberia (figure 1(d)), and Myanmar exhibited no consistent runoff ratio increases (increases observed in 1 of 11 catchments, 9%). The temperate and subtropical climates contained regions with increases (figures 1(b), (c)) and no significant increases (figures 1(f), (g)) (increases observed in 12 of 22 catchments, 54%). Monthly distributions of storm event RR and significance results are presented in figure S2. Phase shift, θ, was similar within each Köppen-Geiger climate region (figure S3). A peaksover-threshold analysis of pre-and post-deforestation precipitation extremes suggested the frequency and magnitude of daily precipitation were largely stationary within each catchment across the study period ( figure S4).

Hierarchical clustering of catchments by climate, tree phenology, area, and forest loss
On each polar plot (figures 2(a)-(d)) θ is represented by the angle, and MAT, min(EVI), ln(Area), and ln(Forest Loss) are represented by the distance from the origin. Hierarchical clustering by θ and MAT closely followed Köppen-Geiger climate classifications (section S6, figure S5) and provided the strongest grouping of similar runoff responses to deforestation by both Purity and NMI (figure 2(a)). Clustering by MAT and median annual minimum EVI provided separation between catchments with no strong annual EVI pattern, though the more complex responses of subtropical and temperate catchments (green) were not well described (figure 2(e)). Temperate and subtropical catchments exhibited significant responses to forest cover loss when annual precipitation extremes lagged the peak growing season (figure 2(a)). Minimum annual EVI highlights a similar seasonality to MAT, but fails to define coherent clusters of runoff responses ( figure 2(b)). Despite previous studies emphasizing regional importance of catchment size and area of loss in predicting hydrologic responses to deforestation, this analysis of 50 catchments with continuous daily streamflow data suggests these are secondary considerations (figures 2(c), (d), (f), (g)).
Hierarchical clustering by θ and MAT (figure 2(a)) lends some interpretation of water yield and land surface energy responses to deforestation (figure 3). The clusters containing predominantly tropical (Am, figure 3(a)), Mediterranean (Csa, figure 3(c)), and cold semi-arid catchments where precipitation extremes lag the growing season (Bsk, figure 3(d)) exhibited decreases in the evaporative index in years following forest loss. In contrast, clusters containing humid subtropical (Cfa, figure 3(b)), humid continental (Dfb, figure 3(e)), and cold semi-arid catchments where precipitation extremes precede the growing season peak (Bsk, figure 3(f)) exhibited a minimal change in water yields following forest cover loss.

Global alignment of EVI and precipitation extremes
MAT is generally aligned by latitude and insolation (figure 4(a)), whereas θ ( figure 4(b)) is more spatially varied, driven by global temperatures, landforms, and regional weather patterns. Spatial patterns of the month of annual maximum precipitation and peak EVI that define θ (figure 4) globally are presented in figure S5. Large-scale mechanisms of extreme precipitation such as the fall season tropical systems along the Atlantic coast of CONUS, winter tropical moisture export precipitation along the Pacific coast of CONUS, Brazilian winter and East African summer monsoonal precipitation, and South African summer precipitation define regional varied patterns of θ across landforms ( figure 4(b)).

Alignment of climate and tree phenology influences runoff responses
Forest cover loss can decrease total catchment transpiration (Farley et al 2005, Filosa et al 2017), and increase catchment saturation during the summer growing-and fall rewetting-seasons (e.g. Knighton et al 2019). Increased saturation (i.e. less opportunity for infiltration and storage of rainfall), will lead to a greater proportion of rainfall becoming runoff. During winter and spring months in subtropical, temperate, and boreal catchments, soils are typically saturated irrespective of the transpiration that occurred during the preceding summer months (e.g. Soulsby et al 2017, Knighton et al 2019). The cold semi-arid (Bsk) catchments of South Africa and Australia showed no consistent increase in runoff ratios to deforestation (figure 2(a)). In these catchments, precipitation peaks occurred approximately three months before peak EVI, during periods in which the magnitude of soil saturation was likely not affected by forest transpiration. In contrast, catchments in the US where precipitation peaks lagged peak EVI (i.e. where forest transpiration would have likely influenced soil saturation at the time of peak rainfall), experienced more consistent increases in runoff ratios (figure 2(a)). Berghuijs et al (2016) proposed that the timing of daily peak discharge across CONUS (characterized by temperate and subtropical climates) is largely determined by trends in catchment storage dynamics and snowmelt. Sivapalan et al (2005) theorized this is because the sensitivity of runoff to storage is greatest when precipitation extremes are aligned with periods where catchment stored water varies (i.e. transitional seasons when soils are not always saturated or dried). Ye et al (2017) modified this by highlighting the interplay between the seasonal alignment of precipitation and energy cycles which together could offer a credible hypothesis of the dominant environmental drivers of daily runoff generation. Our analysis offers an analogous interpretation for the influence of subtropical and temperate forest loss on event runoff: rewetting season runoff is increased in subtropical and temperate catchments where the annual cycle of precipitation extremes lags the annual energy cycle (i.e. precipitation peaks align with tree dormancy and onset of fall rewetting). Importantly, forests can only increase available storage to the physical capacity provided by their canopy and underlying soils. Soil water storage determined by underlying geology likely influences the potential for forest cover loss to increase runoff.
Geographic variations in runoff responses to forest cover change in temperate and subtropical regions (e.g. Preti Trimble et al (1987) examined continuous discharge records of river basins within the humid subtropical Piedmont of the US Atlantic coast following conversion from agriculture to forest. They observed decreases in annual runoff totals following afforestation; however, changes were minimal in wet years when AET losses from saturated agricultural fields were hypothesized to equal that of well-established forests, highlighting the importance of soil water storage. Further in agreement with our findings, Hewlett and Bosch (1984) and Brown et al (2013) found that afforestation of perennial catchments in South Africa yielded changes predominantly to summer season low flows, with limited impact on event-based and multi-year runoff, respectively. Liu et al (2019) suggest decreases in annual runoff coincident with decreased vegetation cover in Southeastern Australia were related to long-term declines in total precipitation, offering an alternative explanation. We observed no consistent change in event runoff ratios in high MAT Monsoonal regions (figures 2(a), S2(a)), possibly owing to a lack of seasonality in insolation, EVI (figure S3), and therefore less annual variation in catchment saturation.
Though our study did not identify catchment size and area of deforestation as significant predictors of changes in runoff (figure 2), the hydrologic importance of these variables has been demonstrated by others , Leite-Filho et al 2019. We propose that future research with larger datasets could consider clustering by additional variables (e.g. MAT, θ, and area of forest loss).

MAT influences runoff responses
We observed consistent increases in spring runoff at low MAT catchments ( figure 2(a)). This may be related to the potential for increased snowpack accumulation following forest cover loss in cold-weather catchments. Increases in snowpack are likely to be followed by increased melt water, soil saturation, and possibly runoff over the subsequent spring. Increased incoming solar radiation (insolation) reaching the snowpack because of canopy loss could have a complicated influence on spring runoff. Increased radiation absorbed by the snowpack could decrease accumulation throughout winter months, reducing spring melt runoff. Conversely, increased insolation could lead to earlier and more rapid periods of melt and greater spring runoff. Our approach focused on rainfall-driven runoff and therefore neglects periods of snowmelt that occurred outside of rain-on-snow events. Soulsby et al (2017) and Geris et al (2015) hypothesized that variations in tree cover would have limited influence on rainfall/runoff partitioning on snow-free soils in northern, low-energy catchments. Our study supports this; however, we note that neither set of authors directly examined the potential for forest cover loss to increase snow accumulation and spring runoff, mechanisms that appear to shift in cold weather catchments (figure S2). Mao and Cherkauer (2009) estimated that deforestation in the humid continental Great Lakes region of the US, led to large increases in spring and summer runoff. They proposed that dense tree canopies provided interception, enhancing sublimation losses, and limiting total ground snowpack accumulation. Loss of canopy shading resulted in greater insolation, more frequent and intense melt events, increased soil wetness, and runoff. Our findings and those of Mao and Cherkauer (2009) suggest that in colder regions with snow accumulation and melt cycles, runoff may increase during periods of snowmelt (figure S2), and are less sensitive to the alignment of precipitation and energy cycles ( figure 2(a)).

Forest loss influence on ET and streamflow partitioning
Many previous studies have indicated the global potential for forest cover loss to increase (and afforestation to decrease) catchment water yield (e.g.  Rice and Emanuel (2019) and Scaife and Band (2017).
Trees influence water partitioning within the critical zone between atmospheric return (i.e. evapotranspiration) and catchment yield (i.e. streamflow). Transpiration by trees can exceed evapotranspiration from nearby non-forested regions if root water uptake (RWU) from deeper soil waters is significant (Brantley et al 2017, Evaristo and McDonnell et al 2017, Barbeta and Peñuelas 2017, Knighton et al 2020. Tree canopies are seasonally dynamic reservoirs that retain precipitation and thus reduce groundwater recharge and surface runoff. In northern low-energy catchments, trees play a critical role in snow accumulation and melt through winter snowfall interception, canopy sublimation, and shielding of the land surface from solar radiation (e.g. Mao and Cherkauer 2009).
Several ecohydrological mechanisms could explain the lack of a forest cover loss effect on catchment EI, exhibited by some climate regions ( figure 3). First, forest cover loss will increase alternative ET pathways. Decreases in EI following a forest cover loss would be reflective of either (1) forest tree water needs exceeding that of successional vegetation or (2) plant-accessible water by forest trees exceeding that of evaporation or successional vegetation. It has been hypothesized that RWU of deeper soil water may be a strategy of certain trees to survive periods of low shallow soil water content (e.g. Matheny  Cavalcante et al (2019) studied catchments in the deforestation arc of the Amazon finding that deforestation events before 1994 resulted in decreased annual EI. After 1994, deforestation had minimal influence on EI, which they attributed to changes in land management practices.
There is also some interplay between catchment storage and water yield responses (Mcdonnell et al 2018). Nijzink et al (2016) showed that deforestation significantly altered rooting zone water storage capacity in a small Mediterranean catchment dominated by Douglas fir, suggesting the influence of trees on catchment water yields includes mechanical modification of soils by roots. Preferential flow pathways may also develop along tree roots, increasing infiltration capacity (Johnson and Lehmann 2006). Within the African drylands, the presence of Vitellaria paradoxa enhanced infiltration for precipitation events exceeding 20 mm 1 day −1 in areas with high soil water storage (150 cm), but had minimal influence on infiltration in regions with lower soil storage (50 cm) (Bargues-Tobella 2020). Peña-Arancibia et al (2019) suggested that the hydrologic response to deforestation was related to climate seasonality, soil water storage, and the regional groundwater recession rate. The storage reservoir provided by the rooting zone must be sufficiently sized to modify land surface runoff responses to forest cover change. In the case of a shallow vadose zone, the sensitivity of runoff to RWU will be physically limited by available rooting zone water storage.

Significance
Global efforts to maintain and increase forest cover require a stronger understanding of hydrologic responses to forest loss. Our findings can inform policy related to the maintenance and conservation of existing forests (e.g. UN-REDD + [United Nations (UN) 2020]), and afforestation (e.g. the Bonn Challenge [Dave et al 2018]). Global climate change is expected on average to increase both precipitation intensity (Fischer and Knutti 2016) and flooding (Trenberth et al 2011). Conservation and afforestation offer the potential to provide a natural sink for atmospheric carbon, slowing global climate change, while potentially also offsetting regional hydrologic land surface responses to more intense precipitation (Ellison et al 2017).
For example, we observe that the estimated regions of land which exhibit favorable phase shifts between EVI and precipitation overlap several of earth's large contiguous (or near contiguous) forests in the tropics, subtropics and temperate zones identified by Hansen et al (2013) including: portions of the Amazon, Gabon, Republic of Congo, and Democratic Republic of Congo, the Appalachians of the Eastern US, forests of the Pacific Northwest (US and Canada), and the Transborder Rainforest Heritage of Borneo (Indonesia) (figure 4). Our findings support the conservation and improvement of these sites as a means of limiting future hydrologic extremes in addition to providing broader ecosystem services. Further, our work broadly demonstrates the importance of boreal forests, which extend across North America, Europe, and Asia.
Several global programs aim to increase forest cover such as the United Nations (UN) Strategic Plan for Forests (UN 2017) and UN-REDD + (UN 2020). These efforts cite flood mitigation as a benefit of increasing forest cover; however, subject-specific policy guidance from the UN is cautious about predicting hydrologic benefits of increased tree planting (FAO 2005, Ellison et al 2017. Guidance generally neglects variations in water uptake potential by different species (e.g. Knighton et al 2019). The high costs and uncertain flood mitigation returns on investment possibly make afforestation and reforestation less attractive mitigation options than assumed. Our research attempts to reduce this uncertainty by furthering our mechanistic understanding of how forests influence runoff generation, guiding policy development.

Acknowledgments
This work was supported by the National Socio-Environmental Synthesis Center (SESYNC) under funding received from the National Science Foundation DBI-1639145.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary information files).