Modeling seasonal vegetation phenology from hydroclimatic drivers for contrasting plant functional groups within drylands of the Southwestern USA

In dryland ecosystems, vegetation within different plant functional groups exhibits distinct seasonal phenologies that are affected by the prevailing hydroclimatic forcing. The seasonal variability of precipitation, atmospheric evaporative demand, and streamflow influences root-zone water availability to plants in water-limited environments. Increasing interannual variations in climate forcing of the local water balance and uncertainty regarding climate change projections have raised the potential for phenological shifts and changes to vegetation dynamics. This poses significant risks to plant functional types across large areas, especially in drylands and within riparian ecosystems. Due to the complex interactions between climate, water availability, and seasonal plant water use, the timing and amplitude of phenological responses to specific hydroclimate forcing cannot be determined a priori, thus limiting efforts to dynamically predict vegetation greenness under future climate change. Here, we analyze two decades (1994–2021) of remote sensing data (soil adjusted vegetation index (SAVI)) as well as contemporaneous hydroclimate data (precipitation, potential evapotranspiration, depth to groundwater, and air temperature), to identify and quantify the key hydroclimatic controls on the timing and amplitude of seasonal greenness. We focus on key phenological events across four different plant functional groups occupying distinct locations and rooting depths in dryland SE Arizona: semi-arid grasses and shrubs, xeric riparian terrace and hydric riparian floodplain trees. We find that key phenological events such as spring and summer greenness peaks in grass and shrubs are strongly driven by contributions from antecedent spring and monsoonal precipitation, respectively. Meanwhile seasonal canopy greenness in floodplain and terrace vegetation showed strong response to groundwater depth as well as antecedent available precipitation (aaP = P − PET) throughout reaches of perennial and intermediate streamflow permanence. The timings of spring green-up and autumn senescence were driven by seasonal changes in air temperature for all plant functional groups. Based on these findings, we develop and test a simple, empirical phenology model, that predicts the timing and amplitude of greenness based on hydroclimate forcing. We demonstrate the feasibility of the model by exploring simple, plausible climate change scenarios, which may inform our understanding of phenological shifts in dryland plant communities and may ultimately improve our predictive capability of investigating and predicting climate-phenology interactions in the future.


Introduction
Terrestrial vegetation responds to local and regional climate and its expression of the water cycle through distinct phenological events, as influenced by the seasonal variations in temperature and precipitation (Cleland et al 2007, Richardson et al 2013. Plants experience the direct and indirect effects of temperature and precipitation through the timing of leaf-on and leaf-off, as well as the seasonal availability of water in the unsaturated and saturated zone respectively, which drives vegetation productivity and composition (Richardson et al 2013, Schlaepfer et al 2017. Phases of green-up, maturity, senescence and dormancy, translate into annual peaks and troughs of greenness, which define the seasonal phenology as an expression of plant growth and vegetative health in response to hydroclimate forcing. Species-level observations have identified temperature as the main driver of spring onset and a break in cold season dormancy in temperate and boreal forests (Hunter and Lechowicz 1992, Hänninen and Kramer 2007, Berra and Gaulton 2021. However the level of climatic controls have been shown to vary between temperate and dry ecosystems Running 2004, Xin et al 2015). Particularly in semi-arid ecosystems, seasonal moisture availability is recognized as the primary controlling resource of vegetation growth and functioning (Rodriguez-Iturbe et al 1999, Rodriguez-Iturbe 2000. The critical link between phenology and climate variability, particularly temperature and precipitation, is considered a powerful indicator of past climate change with observed phenological changes also serving as guidance to predict the causes and consequences of potential future phenological shifts (Menzel et al 2006, White et al 2009, Diez et al 2012, Matthews and Mazer 2016, Munson and Long 2017, Rafferty et al 2020.
To date, prolonged periods of drought and increased water stress have already affected vegetation across a range of biomes, elevational gradients and species: from grasses and shrubs at low elevation (Yao et al 2006, Cui et al 2017, Munson and Long 2017 to boreal forests (de Beurs and Henebry 2010). Furthermore, changes to natural flow regimes and water table elevations, either as a result of temperature and precipitation changes or anthropogenic alterations, have increasingly put pressure on groundwater dependent ecosystems. This has led to increased mortality and widespread vegetation die-off throughout riparian ecosystems (Stromberg and Tiller 1996, Asner et al 2016, Goulden and Bales 2019, Kibler et al 2021, Rohde et al 2021. Similarly, shifts in the seasonal availability of soil moisture and timing of phenological events, as a result of warming temperatures, have affected soil moisture dependent vegetation, leading to prolonged periods of senescence, widespread vegetation die-back and the spread of invasive species (Liu et al 2012, Gremer et al 2015, Wallace et al 2016, Munson and Long 2017, Warter et al 2020. The effects of warming on green-up and senescence, in conjunction with changes to seasonal water availability, have subsequently raised the occurrence of early browning and senescence throughout forest and grassland ecosystems in the past. As a result, the landscape contained an abundance of dried out and easily ignitable fuel, which has amplified the extent and intensity of wildfires in recent decades, a trend which is likely to continue in the future (Westerling et al 2006, Abatzoglou and Williams 2016, Rao et al 2022.
Despite the sensitivity of phenology to climate change and observations of vegetation responses to changing climatic conditions, the magnitude and extent of the effects on vegetation phenology for specific plant functional types remain poorly understood (Cleland et al 2007, Körner and Basler 2010, Richardson et al 2013. Nevertheless, there is a growing urgency to expand and generalize our understanding of temporal vegetation responses to hydroclimate forcing in different plant biomes. This would support a more complete historical understanding of vegetation response to hydroclimate extremes and inherently improve the representation of phenology within ecohydrological and phenological modeling. An improved parameterization of phenology in water balance models would enable a better understanding of future climate-phenology interactions and generally creating a more accurate representation of phenology within climate models and global change studies (Peñuelas et al 2009, Diez et al 2012, Richardson et al 2012, Haynes et al 2019.
Climate change projections regarding the seasonal delivery of precipitation and a warming atmosphere threaten to raise physiological stress levels beyond existing tolerance thresholds, potentially propagating large-scale phenological shifts (Kimball et al 2010, Gremer et al 2015, Walker et al 2015. Given that dryland vegetation dynamics are primarily driven by the availability of root zone soil moisture (Rodriguez-Iturbe et al 1999, D'Odorico et al 2007, D'Odorico and Porporato 2019, plants have adapted by occupying very specific locations in the landscape in response to high air temperature and localized water availability (Gremer et al 2015, Munson and Long 2017, Bradford et al 2020. Rapid response and recovery strategies have been observed throughout dryland grassland ecosystems to maximize favorable moisture conditions after prolonged dry periods (Pennington andCollins 2007, Munson 2013). At the same time, significant changes to natural flow regimes and water table elevations have compromised the sustainability of riparian vegetation communities (Stromberg and Tiller 1996, Stromberg et al 2005, 2013, Rohde et al 2021. The projected trends towards continued groundwater recession and rising evaporative demand in such water-limited environments thus threaten sensitive dryland vegetation. For example, vegetation shifts towards shrubbier, more drought tolerant species, and increased forest mortality, may occur as groundwater access becomes increasingly unavailable to all but the most deeply rooted species (Stromberg and Tiller 1996, Lite and Stromberg 2005a, Stromberg et al 2013, Sabathier et al 2021.
Ultimately, the true extent of these phenological shifts and overall changes to species composition within both groundwater and soil moisture dependent dryland vegetation communities is not yet fully understood. As such, further analysis is required by improving ecohydrological modeling capabilities over a range of spatiotemporal scales (Walker et al 2012, Gremer et al 2015, Munson and Long 2017, Schlaepfer et al 2017. Analyses of vegetation phenology through remote sensing have significantly improved our ability to track contemporary phenology and vegetation cycles on seasonal and interannual scales over a wide range of species and biomes (Choler et al 2010, de Beurs and Henebry 2010, Walker et al 2014, Lu et al 2015, Moon et al 2021. Previous studies have used climate-phenology relationships to build phenology models, however, they frequently focused on meteorological drivers such as temperature, humidity growing-degree days or vapor pressure deficit (Farrar et al 1994, Ji and Peters 2004, Jolly et al 2005, Choler et al 2010. Few also consider the whole seasonal cycle, but instead focus on individual phenological transition dates, such as leaf-out or senescence (Mariën et al 2019, Chen et al 2020. This is insufficient, particularly in dryland ecosystems, where seasonal water availability instead of temperature is limiting phenology (Richardson et al 2013, Renwick et al 2019. As such, there remains a need for models that capture the temporally complex vegetation dynamics of more water-controlled environments in a comprehensive way. Therefore, in this study we use remotely sensed vegetation indices together with contemporary hydroclimate data to resolve the linkages between climate variability and phenology on a seasonal scale. The main goal of this study is to parameterize interannual variations in greenness, to support ecohydrological modeling within simple water balance models, which typically use fixed annual greenness values that are unresponsive to changes in climate and water availability. Such models often do not even incorporate the seasonality of deciduous vegetation in terms of its influence on the water balance (ET losses via transpiration should not be occurring when the plants are dormant-after leaf-off). Therefore, our modeling goes a step further because it not only captures the effects of green and senesced seasonal periods, but it allows values of greenness to respond to hydroclimatic forcing, such that, for example, high rainfall delivery in the right period will lead to a strong impact on greenness and vice versa.
Overall, we focus on different plant functional types occupying distinct rooting locations within a dryland ecosystem in SE Arizona, from riparian trees to semi-arid grasses and shrubs. The overall objectives of this study are to (a) identify seasonal phenology patterns of different plant functional groups by exploring historical climate-phenology relationships; (b) determine key phenological events along the seasonal phenology curve; (c) develop empirical models for each key event between hydroclimate variables and vegetation greenness and (d) create synthetic seasonal phenology curves that capture vegetation greenness in response to hydroclimate forcing. We will delineate the relative controls of water availability and evaporative demand, including temperature, on phenology of four plant functional groups and use our model to explore vegetation responses under simple, plausible climate change scenarios. Although the results are specific to SE Arizona, our intent is to present a modeling framework that can be applied and tested in other dryland environments to evaluate hydroclimate-phenology relationships and the potential effects of climate change on vegetation dynamics.

Study region
We used remote sensing to sample vegetation greenness of different plant functional groups within two distinct and complementary areas in SE Arizona; the Walnut Gulch Experimental Watershed (WGEW) and San Pedro Riparian National Conservation Area (SPRNCA). The San Pedro and Walnut Gulch are both well studied dryland vegetation environments due to their characteristic hydrology and vegetation dynamics (Renard et al 2008, Skirvin et al 2008, Singer and Michaelides 2017, Stromberg et al 2017, Mayes et al 2020, Sabathier et al 2021. The hydrology of SE Arizona is strongly driven by the North American monsoon. Precipitation patterns exhibit a characteristic bimodal distribution with dual precipitation peaks in winter/spring and summer (Nichols et al 2002). Annual precipitation totals fall in the typical range between 300-400 mm year −1 , with two-thirds of the annual precipitation and 90% of its runoff delivered during the summer months between July and September (Nichols et al 2002, Stromberg et al 2006. Riparian vegetation distribution and dynamics along the San Pedro River are strongly shaped by the spatial variability of available moisture, topography, as well as underlying geology and elevation (Makings 2005, Stromberg et al 2005, Sabathier et al 2021. Based on hydrology and geomorphological characteristics, the San Pedro has been divided into reach types that summarize flow conditions and associated vegetation communities along the river (Stromberg et al 2006). The major vegetation communities along the San Pedro include riparian woodlands on low floodplains, xero-riparian trees and shrubs on river terraces, and grasses and desert shrubs on piedmont surfaces further upland (Makings 2005, Sabathier et al 2021. The narrow gallery forest along the river channel is primarily populated by obligate phreatophytes such as Fremont cottonwood (Populus fremontii) and Goodding's willow (Salix gooddingii), which prefer consistent access to groundwater and are therefore strongly influenced by variations in connectivity between surface and subsurface water flows (Makings 2005, Stromberg et al 2017. Bordering the riparian corridor are upland terraces dominated by mesquite woodlands (Prosopis velutina), which have been observed to be less dependent on groundwater and more flexible in their use of seasonally available water sources, including soil moisture (Snyder andWilliams 2000, Stromberg et al 2017).
We sampled the vegetation of the phreatophyte-dominated riparian floodplain and mesquite woodland terrace longitudinally along a 28 km stretch of the river ( figure 1(a)). This includes wet, intermediate and dry riparian flow reaches, spanning a gradient from perennial to intermittent flow (Stromberg et al 2006, 2017, Sabathier et al 2021. The delineation of reaches is based on stream flow hydrology and geomorphological properties. It is also based on vegetation bioindicators that are sensitive to surface water or groundwater and has been established by Stromberg et al (2006) and is widely used to differentiate vegetation and hydrological conditions throughout the SPRNCA.

Hydroclimate
We downloaded several datasets of climate variables to support our analysis of hydroclimate controls on phenology. Daily precipitation and minimum/maximum air temperature data for the SPRNCA were extracted from the National Oceanic and Atmospheric Administrations (NOAA) Climate Prediction Center's (CPC) Unified Gauge-Based Analysis of Daily Precipitation (https://psl.noaa.gov/data/gridded/data.cpc. globalprecip.html) and Global Unified Temperature (www.psl.noaa.gov/data/gridded/data.cpc.globaltemp. html). The CPC data combines the available gauge-based precipitation information and uses an optimal interpolation objective analysis technique to create a gauge-based spatial precipitation output (M. Xie et al 2007, Chen et al 2008. We generated datasets for the period 1994-2019 with a spatial resolution of 0.5 • × 0.5 • covering our study area along the San Pedro River. At Walnut Gulch, precipitation data were extracted from the Southwest Research Center Data Access Project (www.tucson.ars.ag.gov/dap/) for the period between 1994-2019. Precipitation data is available for a multitude of gauges located throughout the catchment. However, we used only a selection of the available rain gauges (Nr: 1,3,4,5,7,8,45,46,48,50,53,56,57,58,59,62,63,64,67,72), which are closest to the vegetation sampling points of grass and shrub (figure 1(a)). We then created a spatial arithmetic mean between all gauges for our analysis. It should be noted that there was no significant difference between the stations, as they were all relatively close to each other. At Walnut Gulch we were able to get local precipitation data from the well-studied USDA network of rain gauges in the catchment. However, at the San Pedro Riparian Area, local climate data is not readily available, thus the NOAA climate dataset at 0.5 • was chosen as the most suitable. As such, we do not expect dramatic differences in the seasonal distribution of rainfall in this system, and the strong hydroclimatic correlations between vegetative greenness and rainfall bear this out for the different plant functional groups. Ultimately, our goal was to create a general model structure that could be applied in a variety of settings, which might also face similar data availability challenges.
Daily potential evapotranspiration was extracted from the hPET dataset for our two study sites (https:// data.bris.acuk/data/dataset/qb8ujazzda0s2aykkv0oq0ctp). To characterize water availability to phreatophytes we obtained datasets of mean daily depth to groundwater (DTG) from the US Geological Survey's National Water Information System (https://waterdata.usgs.gov/az/nwis/), focusing on monitoring wells with the most continuous measurements and longest monitoring periods. We obtained records of several groundwater wells (USGS ID: 313 738 110 102 901 (GS_SP17), 314 511 110 120 601 (GS_SP23) and 314 904 110 125 001 (GS_SP27)) located in the wet, intermediate and dry reaches respectively (figure 1(a)). DTG records for wells GS_SP17 and GS_SP23 contain daily data between 2001-2018, with 80% and 52% of the respective daily records complete. This means that data gaps exist in the timeseries, with a certain number of days missing. Therefore, a fewer number of data points could be used to assess DTG at the time of maximum greenness for riparian terrace and floodplain trees. GS_SP27, being a periodic monitoring well, only has limited observations available between 2007-2010. Therefore, we did not use this well in our analysis due to the limited amount of data available. DTG generally varies among reaches due to local differences in the underlying geology and geomorphology, as well as proximity to tributaries and ground-water pumping sites (Lite and Stromberg 2005a). In this study, we primarily explored hydroclimate-phenology interactions for the perennial and intermittent reaches. We chose to frame our analysis within the range of a water year, which generally refers to a time period of 12 months for which precipitation totals are measured. It is usually timed to begin after a dry period, so the moisture input can be 'reset' . In many hydrological studies, it is set to start on October 1st, with days counted sequentially from the first day of the water year (DOWY) (i.e. October 1 st ) until the end (i.e. September 30 th ). In this study we set the water year from November 1 st to October 31 st to represent the bimodal precipitation regime of this monsoon dominated ecosystem and the associated vegetation dynamics of the region.

Satellite remote sensing of vegetation greenness
In order to measure how vegetation greenness evolves seasonally and responds to inter-annual variations in climate, we obtained continuous timeseries of surface reflectance images from the Landsat 5 Thematic Mapper, and Landsat 8 Operational Land Imager/Thermal Infrared Sensor (https://espa.cr.usgs.gov/). The Soil Adjusted Vegetation Index (SAVI) is considered to be a particularly suitable metric in regions with relatively low vegetation cover (Huete 1988). We only selected cloud-free images which are available at 30 m spatial resolution every 16 d. SAVI is in general closely related to the Normalized Differenced Vegetation Index (NDVI) but considered less sensitive to soil brightness. The effects of soil background reflectance in areas where vegetation cover is less dense are accounted for by including a soil brightness correction factor (L) (Huete 1988). SAVI is defined as: where L is the soil brightness correction factor, commonly set to 0.5, which is considered applicable to suit most land cover types (Huete 1988).
To avoid mixed SAVI signatures of multiple vegetation types at Walnut Gulch, we visually identified homogenous areas where grasses and shrubs are the dominant vegetation cover type by using a map of vegetation classes previously established by Skirvin et al (2008), as well as most recent aerial imagery from the National Agricultural Imagery Program (NAIP) at a 1 m spatial resolution (www.fsa.usda.gov/programsand-services/aerial-photography/imagery-programs/naip-imagery/index). To validate the selected SAVI pixels we checked for potential changes in land cover/ vegetation type using maps from the USGS National Land Cover Database, which is updated every five years, for the period 2001-2019 (www.usgs.gov/centers/ eros/science/national-land-cover-database). If no significant changes in the land cover change index was found, it can be assumed that the land cover type did not significantly change during the study period. For Landsat images, we used a single pixel approach, selecting 30 individual sampling points clustered (cloud) over a homogenous area in the Walnut Gulch Catchment. For each point we extracted remote sensing data (SAVI) between 1994-2021. Similarly, at the San Pedro site we used high-resolution NAIP imagery to visually delineate the riparian floodplain and mesquite terraces and sampled multiple pixels for each plant functional type longitudinally along a 28 km stretch ( figure 1(a)). To obtain the most homogenous signature of vegetation greenness for the riparian floodplain, we selected only pixels in areas where the width of the cottonwood corridor exceeded 30 m.
We further filtered the obtained SAVI data to include only years where images in June, July and August were available, which usually corresponds to the period when riparian vegetation is leafed out, and plants exhibit maximum greenness values. The SAVI from all sampled pixels was summarized into a spatial median for each plant functional group over the study period. We then aggregated the SAVI with its 16 d frequency to mean monthly values to generate 12 monthly averages per year. Each month was then averaged across years (n = 27) into a seasonal annual composite (i.e. every January over the study period). This leaves a composite series of 12 values for every month, including standard deviation of interannual variations. This not only smooths the SAVI signal, but also removes any issues regarding aligning sampling dates from remote sensing data with climate data. We use the seasonal composite to characterize the phenology curve and determine the key phenological events and controls on timing and amplitude of greenness.

Phenological model development
The framework and step-wise approach for the proposed methodology of creating synthetic phenology curves is broadly illustrated in figure 2. To model the key controls of seasonal greenness of different plant functional groups, we first analyzed the statistical relationships between hydroclimate and remote sensing data. We first defined the key events along the seasonal curve and then extracted the greenness during each key event for every year. We then explored the possible link to relevant hydroclimate variables, by establishing empirical relationships between greenness at the time of the key event (i.e. maximum greenness) and a relevant hydroclimate variable. We did not create long-term mean seasonal climate variables but used the annual climate data for the regression analyses. Regression models were developed based on relationships between hydroclimate and remote sensing data and used to predict the amplitude of greenness during each event.

Key phenological events
For defining the key phenological events for all plant functional groups we considered the shape of the mean seasonal SAVI composites. The seasonal cycle of vegetation in SE Arizona can show a characteristic unimodal or bimodal cycle in response to the local hydroclimate and water availability, and depending on the plant  functional type (Jenerette et al 2010, Notaro et al 2010. Grasses and shrubs generally displays a bimodal seasonal cycle which can be attributed to seasonal interactions with temperature and precipitation. The plants' sensitivity to soil moisture drives the dual greenness peaks, with a first small seasonal peak in spring (POS_1) in response to a break from cold season dormancy mediated by temperature and the accumulation of cool season precipitation (Notaro et al 2010). The increased transpiration and declining moisture resources incur a period of low photosynthetic activity before the arrival of the monsoon. This leads to a second larger seasonal peak (POS_2), after which greenness starts to subsequently decline on a steady downwards trend until the end of the season (EOS) (figure 3(a)). In contrast, the unimodal cycle of riparian trees exhibits a monotonic green-up in spring and a single seasonal peak at the height of the growing season between July and August (figure 3(b)) (Notaro et al 2010, Walker et al 2014. We defined several key phenological events based on the two characteristic phenology distributions (unimodal and bimodal), which we use as hinge points to reconstruct the seasonal phenology curves. The start of green-up (SGU) for all plant functional groups is defined as the DOWY on which the median SAVI shows consistent increase and upwards trend, representing a break in dormancy and the onset of photosynthetic activity. For bimodal plant functional types (figure 3(a)) we define a spring greenness peak (POS_1). The next hinge point is the start of the monsoon (SOM) from where a steady increase in greenness leads to the second larger greenness peak at the height of the growing season (POS_2). After this point, plants can be considered to start senescing as greenness starts to steadily decline until the EOS, which is estimated as the last DOWY (DOWY 365). In unimodal vegetation, after green up they only exhibit a single seasonal peak (POS_2), after which senescence sets in and greenness declines towards the end of season. For certain trees, the peak of season (POS) shows a broader peak, indicating a more extended period of maximum greenness and maturity. In this case, an additional hinge point is defined to characterize the start of the peak greenness period (ps). . Derivatives of mean maximum daily air temperature (red) and mean seasonal composite SAVI (grey) for (a) riparian trees, (b) grasses and shrubs. The zero crossing points of SAVI and temperature indicate the slope changes from negative to positive (and vice-versa) and determine changes in greenness and seasonal temperature respectively. The dates of the zero crossing points determine the timing of SGUt and POSt. ± 1STD is shown as red and grey shading around the composites of air temperature and SAVI respectively.

Timing of key phenological events
The first step in creating a synthetic phenology is to establish the timing of key phenological events based on greenness responses to hydroclimatic drivers. Once these are determined and greenness values are calculated for them, a synthetic curve can be fit via splines ( figure 3). Previous studies observed a close link between the intra-seasonal timing of phenological events and temperature, as well as precipitation. However, temperature is typically considered the main driver of the timing of green-up and senescence (Jolly et al 2005, Richardson et al 2013, Renwick et al 2019. Therefore, we have adopted that concept here to define the start of the phenological green up (SGU t ) and to indicate the timing of the highest peak of season greenness (POS_2 t ). To extract the phenological timing for these two events, we utilize SAVI and temperature data. Given that SAVI is only measured every ∼16 d, but temperature is measured daily, we had to account for the time lags between temperature and SAVI by applying the delayed moving average method (Wu et al 2021) to the timeseries of daily maximum temperature, by using 16 d averages to match the temporal resolution of SAVI. Thus, we developed synchronized curves of both SAVI and temperature from which to analyze the timing of phenological metrics (figure 4).
Next, to identify the timing of SGU and POS for all plant functional groups as a function of temperature, we computed first derivatives of the time-matched curves of SAVI and temperature via the curve derivative method (Zhang et al 2004). Then we quantitatively compared seasonal SAVI derivatives with the derivative of temperature to identify when the slopes of both curves cross zero (change of sign, for example, where a change from negative to positive values indicates an increase in greenness associated with the SGU). The timing of the peak of season (POS_2) is determined as the date when the SAVI derivative crosses zero again (from positive to negative), indicating the start of a monotonic decline in greenness, defined as the start of senescence (figure 4), which continues until the EOS and which is defined as the end of the water year set to October 31 st . Note the relationships between these metrics in derivative space and how they compare to the actual SAVI curves (c.f., figures 3 and 4).
These steps can be carried out to determine the timing of phenological events to model a unimodal synthetic phenology curve, for example that of riparian terrace trees such as mesquite. However, some plant species exhibit more complex phenologies that require additional information. For example, riparian floodplain trees such as cottonwoods in the Southwest USA may have a prolonged, broader peak greenness that needs to be delineated by an additional hinge point, which denotes the start of the period of maximum greenness. The timing of this event is also forced by temperature, as the date when the temperature derivative exhibits maximum change (see figure 4).
For bimodal plant species with two seasonal peaks in greenness, such as for grasses and shrubs, we separately compute the timing for each peak (POS_1 and POS_2, figure 3(a)). However, for these plant functional groups, we found that the temperature derivative curve does not help in defining the timing of the first spring peak, which is associated more with antecedent precipitation, which defines the store of root-zone water availability in these water-limited environments (Jolly et al 2005, Tang et al 2015. Because the bimodal phenology is closely linked to the dual precipitation distribution of the region, we determined the timing of the first peak, POS_1, and the SOM, using the derivatives of seasonal precipitation instead of temperature. SOM t is determined by the onset of the monsoon, as the inflection point with the steepest increase along the cumulative seasonal precipitation composite mean ( Figure S1(a)). POS_1 t is determined using the derivative of precipitation at the point where the derivative equals zero, which indicates the end of the spring rains and thus the peak of the first period of green-up ( Figure S1(b)).

Greenness amplitude for key phenological events
To estimate the sensitivity of different hydroclimate variables on greenness amplitude, we performed correlation analyses between hydroclimate and observed SAVI values of key phenological events. We extracted annual values for each phenological event and tested the correlation strength between observed SAVI and concurrent and antecedent/lagged hydroclimate variables. Given the direct sensitivity of soil moisture dependent species, such as grasses and shrubs, to accumulated antecedent precipitation in the soil (Laio et al 2001, Porporato et al 2002, D'Odorico et al 2007, we applied the moving average method to the timeseries of daily precipitation. Using multiple 16 d composite periods we tested lags of up to 96 d to account for antecedent rainfall contributions to soil water availability (Ukkola et al 2021).
Vegetation dynamics of deeper rooted phreatophyte species such as riparian floodplain trees, are known to respond primarily to variations in water table depth, depending on root zone access to the water table (Porporato et al 2002, Laio et al 2009, Stromberg et al 2017. However, it has also been shown that riparian trees may access water from the upper soil layers to substitute water demands (Snyder andWilliams 2000, Singer et al 2014). This soil water availability might be reasonably generalized as antecedent available precipitation (aaP = P − PET), which we use a proxy for deeper root zone water availability. Therefore, we analyzed both depth to groundwater and aaP as potential controls on the seasonal greenness peak for riparian corridor and terrace trees, testing the relationships between SAVI and DTG and between SAVI and aaP. Due to the unimodal phenology cycle in riparian trees, we determined POS as the most responsive phenological event to hydroclimate variations and created normal distributions based on observed values for greenness at SGU and EOS, as these showed minimal interannual variations.
To account for the differences in the DTG of different wells, observed data were normalized to the same range and scale (0-1), irrespective of absolute values. DTG values are expressed as DTG_norm and each well is normalized as: Based on the functional form of regressions between observed SAVI and the relevant hydroclimate variable for each phenological event, we fit linear regression models, which were calibrated using half of the available data. We used antecedent precipitation (P) for grasses and shrubs and depth to water table as well as antecedent available precipitation (P-PET) for riparian floodplain and terrace trees. Specifically, we calibrated the regression models to a subset of data using odd or even years between 1994-2021 (n = 13 years) and assessed the goodness of fit using either Pearson's R or Spearman's Rho correlation. For the purposes of creating the most dynamic synthetic phenology curves, we took into account the inherent uncertainty in these relationships, by introducing a prediction interval (PI) around all regression models. The 95% PI denotes 95% confidence that this area will contain a new observation of greenness in response to a climate forcing and is defined as: where y pred is the predicted value of greenness in response to climate forcing x * , T denotes the 97.5th percentile of the Student's t-distribution with n-m degrees of freedom, and σ is the standard deviation of the residuals at the confidence level (Bevington and Robinson 2003). We then used the prediction interval as a range for a Monte Carlo sampling loop to create multiple simulations (N = 100) of stochastic greenness values that would represent the spatial and temporal variability of observed greenness values in response to hydroclimate forcing at the temporal location of each phenological event.

Modeling of synthetic phenology curves
To create continuous synthetic phenology curves, we used the coefficients from the regression models of each phenological event to estimate annual synthetic SAVI values at the timing of a key phenological event in response to hydroclimate forcing. Stochastic simulations of greenness (N = 100) for each event within the range of the prediction interval are summarized into a sample median. We validated the performance of the regression models and the models' ability to capture greenness responses to variable hydroclimate forcing, by using the other half of the hydroclimate data, comprising the remaining years between 1994-2021 not used for calibration (n = 13 years). We used monotonic cubic spline interpolation between phenological events of an individual year and fit piecewise cubic polynomials, creating smooth continuous functions which pass through all data points without inducing unwanted oscillations in the fitted curve and insure monotonicity (Wolberg and Alfy 2002). The cubic spline are fitted at an interval between phenological events, with the number of data points and intervals varying between bimodal (n = 5) and unimodal (n = 3) phenology cycles (figure S2). We use the PchipInterpolator function in python (https://docs.scipy.org/doc/scipy/ reference/generated/scipy.interpolate.PchipInterpolator.htm) to automate the process of fitting splines between events for every year. The interpolant uses 1-D arrays of x and y, representing time and synthetic SAVI respectively, as well as optional parameters extrapolate, which we set to True. The synthetic SAVI timeseries over the validation period were summarized into a seasonal mean composite, which was compared to the observed SAVI composite of each plant functional group. The correlation between observed and synthetic composites was assessed for using Spearman's Rho correlation while overall distributions were compared using the two-sample Kolmogorov-Smirnov statistic at a significance level of α = 0.05. Additionally, to determine whether a phenological model was able to capture the variability of observed SAVI, we assessed synthetic greenness values of each key phenological event through Kruskal-Wallis testing at a significance level of α = 0.05. We also evaluated timeseries of mean synthetic and observed values through Nash-Sutcliffe efficiency (NSE), as well as the coefficient of determination (R 2 ) to estimate the measure of the variance explained by the model. Overall model performance was deemed adequate if no significant statistical difference between modeled and observed values was detected.
To further test and demonstrate the functionality of the model and explore the potential impacts of climate change on phenology, we created simple plausible climate change scenarios that include changes to key hydroclimate variables. Wintertime precipitation is projected to decline, while monsoon precipitation projections remain uncertain, with some regional climate models pointing towards a drying in monsoon regions of the Southwest, in addition to an overall drying trend and increased drought risk (Cook et al 2015, Pascale et al 2017, McKinnon et al 2021. Based on these projections, we applied a simplified scenario of reduced precipitation intensity by 20% during winter/spring (January-April) as well as during the monsoon (June-September) without changing the frequency of events. We further applied a reduction of PET by 20% to represent the rise in evaporative demand. As groundwater levels have been observed to decline in riparian areas throughout the Southwest (Stromberg and Tiller 1996, Kibler et al 2021, Rohde et al 2021, we implemented a lowered water table depth by 0.5 m to explore differential vegetation responses to a deeper water table. We chose to look at a lowered DTG independently of climate, as a representative of the potential effects of lowered precipitation and increased PET. Furthermore, we also explored the effects of temperature shifts and earlier warming, a trend which has been observed to affect vegetation phenology and advance the onset of green-up as well as senescence (Xin et al 2015, Munson and Long 2017, Warter et al 2020. We examined the impacts of such changes on vegetation responses of all plant functional groups, comparing mean growing season greenness distributions (March-September) between historic observations and model results through two-sample Kolmogorov-Smirnov testing and discuss the potential implications on vegetation and ecosystem functioning.

Timing of phenological events
A comparison of derivatives of composite SAVI and composite maximum air temperature predicted the timing the start of green-up (SGU t ) in riparian floodplain trees in early February (∼DOWY 110), while in terrace trees it occurred in early March (∼DOWY 125) (figures 4(a) and (b)). The timing of the peak of season (POS t ) was similar for trees, occurring around the same time in late August (∼DOWY 300), linking the timing of SGU and POS to air temperature. Grasses and shrubs showed similar links between air temperature and SAVI regarding timing of SGU and the POS. SGU t in spring occurred in early March (∼DOWY 124), while the main summer peak POS t occurred at the end of August (∼DOWY 300), correlating well to the seasonal changes in maximum air temperature. The timing of SOM generally occurred between mid-June to early July (∼DOWY 245) while POS_1 t in shrubs and grasses usually occurred post spring rains and prior to the onset of the early dry season around mid-April (∼DOWY 170).

Greenness responses to hydroclimate forcing
We found significant correlations between DTG and maximum greenness of riparian floodplain (p < 0.05) and terrace trees (p < 0.05) as well as antecedent available precipitation (p = 0.012 for floodplain; p = 0.015 for terrace trees) (figure 5). The groundwater well in the perennial reach is in close proximity to the river (<20 m), with a shallow annual mean DTG of 1.8 m and small inter-annual fluctuations (±0.2 m). The significant negative linear correlation observed between maximum SAVI and DTG, could be well explained through the regression models (R 2 = 0.71, RMSE = 0.01 for floodplain trees; R 2 = 0.72, RMSE = 0.02 for terrace trees) (figures 5(e) and (f)). At the same time, the significant positive linear trend between aaP and maximum SAVI, could also be well explained by a regression model (R 2 = 0.91 for floodplain, RMSE = 0.1; R 2 = 0.82,RMSE = 0.13 for terrace trees) (figures 5(g) and (h)). Overall, mean POS values were similar in both riparian floodplain and terrace trees, however we observed a slight upwards trend over the observation period, despite short term dry periods between 2018-2020 (figures 5(a) and (b)).
Overall, the presence of significant relationships to two different water sources makes sense, as riparian forest sites have been observed to also use shallow soil water during the rainy season (Snyder and Williams 2000). At the intermediate reach, we found the same significant relationships between POS and DTG, which could also be reasonably explained through the regression models (R 2 = 0.75, RMSE = 0.08; R 2 = 0.81, RMSE = 0.003) as well as aaP (R 2 = 0.72, RMSE = 0.03 for floodplain; R 2 = 0.73, RMSE = 0.08 for terrace trees) (figures S2(a) and (b)). The water table was moderately shallow with a mean annual depth of 2.8 m and larger interannual fluctuations (±0.9 m). We ultimately used the strong linear relationship between DTG and POS to establish synthetic phenology curves, as DTG is the more direct measure of vegetation greenness in groundwater-dependent riparian ecosystems. However, where groundwater data are unavailable, aaP can still be used as an alternative to force greenness responses.
In grasses and shrubs, POS showed visibly larger interannual variations than the POS in riparian trees, due the strong link to recent precipitation and the dual precipitation cycle ( figure 6). Mean values of POS were similar in grasses and shrubs (0.22 and 0.20, ±0.05 respectively) and showed a strong variability in response to interannual precipitation differences (figures 6(a) and (b)). More specifically, lower greenness values can be seen in years of lower antecedent moisture and weaker monsoon (<50 mm), such as in 2003 and 2004, during which >50% of Arizona was under extreme drought, as well as during summer 2020, which was an exceptional drought year throughout the Southwestern US (U.S. Drought Monitor, https:// droughtmonitor.unl.edu/DmData/TimeSeries.aspx). The highest POS was observed in response to larger monsoon totals such as in 1999 and 2000 (>140 mm) as well as in 2017, which had a particularly strong monsoon (>200 mm). Overall, POS was significantly related to antecedent monsoon precipitation of the previous month (p < 0.01 for grass and shrubs). The interannual greenness variations during the POS could be well explained through the regression models (R 2 = 0.86, RMSE = 0.07 for grasses; R 2 = 0.88 shrubs, RMSE = 0.05, figures 6(e) and (f)). Relationship to DTG for both, grasses and shrubs were not evaluated, because no DTG data was available for this area. Antecedent precipitation was the most significant variable also for other key phenological events. POS_1 was significantly related to accumulated moisture from weak synoptic rain events during the cool season (p < 0.01) (figures S4(a) and (b)), while SOM correlated to the limited precipitation during the early summer dry period (p < 0.01), during which time only limited increases in greenness occurred (figures S4(c) and (d)). Finally, greenness at EOS was most significantly related to antecedent precipitation from late monsoon and early autumn precipitation (p < 0.01) (figures S4(e) and (f)).

Modeling synthetic phenology curves
Synthetic greenness values at POS forced by DTG for riparian floodplain and terrace trees showed no statistically significant difference to observed values (p = 0.81, stat = 0.3; p = 0.62, stat = 0.25), as did values at SGU (p = 0.07, stat = 0.24; p = 0.83, stat = 0.06) (figures 7(a) and (b)). The model showed strong predictive abilities, with a NSE of 0.94 for riparian floodplain and terrace trees between synthetic and observed timeseries (figures 7(c) and (d)). Median composite values showed strong correlation with In grasses and shrubs, synthetic values of POS_1 forced by antecedent precipitation showed no statistically significant difference to observed SAVI (p = 0.6, stat = 0.28; p = 0.89, stat = 0.02). The timing of the SOM was accurately captured and greenness values at SOM showed no significant differences to observations (p = 0.93, stat = 0.006; p = 0.81, stat = 0.06) as well as similar interannual variability in response to variable precipitation totals during the early summer dry period (figures 8(a) and (b)). Synthetic POS values showed interannual variations of greenness in response to variable monsoon totals within the historic range and annual POS values showed no significant differences to observed values (p = 0.62, stat = 0.25; p = 0.9, stat = 0.01). Model performance was confirmed with synthetic and observed composite medians showing no statistically significant differences for grasses and shrubs, as well as a high coefficient of determination (R = 0.81 grasses; R = 0.77 shrubs). Good predictive power was confirmed through high NSE coefficients (NSE = 0.9 grass; NSE = 0.86 shrubs) (figures 8(c) and (d)). Table S2 contains a summary of goodness-of-fit statistics for all key phenological events and plant functional groups.

Vegetation responses to plausible future climate change scenarios
In riparian trees a drop in DTG of 0.5 m led to a significant reduction in annual POS values of approximately 17% in riparian floodplain and 13% in riparian terrace trees (p = 0.02, stat = 0.57; p = 0.02, stat = 0.6, respectively) (figures 9(c) and (d)). Composite means further highlight the effects of shifted temperature, whereby an early onset of SGU and POS is visible. The shift towards an earlier POS, acts as a two-way effect on seasonal greenness as DTG may still be at its seasonal low level, prior to the onset of the monsoon and increase in baseflow. Growing season distributions show significant differences for riparian floodplain and terrace trees (p = < 0.05, stat = 0.37; p = < 0.05, stat = 0.32, respectively), highlighting the shift towards lower POS values (figures 9(c) and (d), inset).
In grasses and shrubs, a 20% reduction in monsoon precipitation was sufficiently large to result in a significant decrease in greenness of POS by 15% in grass and 13% in shrubs (p = 0.03, stat = 0.67; p = 0.03, stat = 0.66, respectively). Similarly, the reduced spring precipitation resulted in more muted green-up and reduced POS_1 (figures 10(c) and (d)). The seasonal composite mean highlights the effects on seasonal phenology patterns, showing the loss of a clear spring green-up peak and shift towards a more prolonged   period of senescence during spring, while the shift in temperature resulted in an earlier timing of SGU and POS (figures 10(e) and (f)). The change in growing season greenness between March and September was strongly driven by reduced antecedent precipitation (R 2 = 0.77 for grass, R 2 = 0.81 for shrubs) and a significant difference between historic and observed (p = < 0.05, stat = 0.18 for grass;p = <0.05, stat = 0.17 for shrubs), with a clear trend towards lower POS and an increased frequency of lower greenness values in spring and summer (figures 10(e) and (f) inset).

Discussion
In this study, we demonstrate the feasibility of using remote sensing and hydroclimate data for the development of a simple, empirical phenology modeling approach for different dryland plant functional groups. In the context of shifting spring green-up and flowering times due warming temperatures, there is a need to prioritize high-resolution temporal and spatial analyses of climate drivers on plant phenology to improve our understanding of climate-phenology interactions that might be expected in the future (Zhang et al 2003, Ganguly et al 2010, Moon et al 2021. The mechanisms of bimodal seasonal phenology cycles have previously been explored in regions where the North American monsoon dominates, however there is remaining uncertainty associated with the complex ecohydrological relationships in drylands of the Southwest USA.
Our results have illustrated the relationships between vegetation greenness and different hydroclimate drivers (figure 11(a)) and their functionality within an empirical modeling framework to explore synthetic greenness responses under plausible climate change scenarios. Such modeling capabilities open up avenues for exploring climate-phenology feedbacks under future climate change scenarios and support the proactive development of appropriate adaptation strategies to mitigate the effects of warming temperatures, increased precipitation variability, and declining groundwater (Joyce et al 2013.
The strong dependence of riparian floodplain and terrace vegetation on groundwater is well appreciated and has been extensively explored with renewed urgency throughout the Southwest and California, due to rapidly declining groundwater levels and increased forest mortality (Stromberg et al 2006, Pettit and Froend 2018, Rohde et al 2021, Sabathier et al 2021. While individual drought years (i.e. 2020) had a relatively strong effect on grasses and shrubs (figures 6(a) and (b)), riparian trees remained unaffected. We attribute this trend to the relatively stable groundwater level and minimal seasonal fluctuations, providing trees with relatively constant access to groundwater to satisfy their transpiration demands despite seasonally limited precipitation and soil moisture. Due to the lack of groundwater observations during the dry years for this well location, we cannot verify this interpretation, however previous studies have shown that flow in the perennial reach of the San Pedro river is consistently maintained by the presence of shallow bedrock and a strong connection between the riparian aquifer and the river, thus buffering the effects of short-term dry periods (Stromberg et al 2006, Baillie et al 2007, Sabathier et al 2021. Ultimately, the survival of riparian ecosystems in the future will depend entirely on the plants' ability to continue deep root growth to match declining water tables, as they are unlikely to survive solely on root zone moisture should groundwater become unavailable (Snyder andWilliams 2000, Pettit andFroend 2018). Considering, for example, the effects of a climatic shift that would lead to a deeper water table and temporal shifts in the timing of green-up, on greenness responses in riparian floodplain and terrace trees ( figure 11(b)); through the critical link of vegetation greenness to groundwater table depth, we found that under a scenario of a deeper water table and earlier green-up, we might expect significant changes to annual greenness values and the timing of phenological events. Subsequently, any associated changes to groundwater-flow patterns as a result of increased precipitation variability, may propagate into reduced stream flow patterns within perennial reaches and risk a shift from perennial to more intermittent flow conditions. Such changes to the seasonal availability of water would have far reaching implications for the health and survival of riparian ecosystems. As such, shifts from dense riparian vegetation towards shrubbier more drought tolerant species might occur, as groundwater becomes increasingly inaccessible or unreliable, ultimately compromising a range of ecosystem functions and risking an increased loss of sensitive habitat and species diversity.
Similarly, shallow rooted species are vulnerable to the potential impacts of precipitation changes and prolonged periods of drought, due to the strong link between antecedent precipitation and seasonal phenology ( figure 11(a) phenological responses into individual events and quantifying the seasonal and interannual variations, we have revealed further details about the precipitation-phenology relationships for grass and shrub plant functional groups and explored their potential vulnerability to future climate change. Simulations of simple climate change scenarios showed that a reduction of spring and monsoon totals could inherently alter phenological patterns and seasonal greenness responses (figure 10). Reduced spring precipitation pulses would lead to more muted spring green-up due to reduced moisture availability, while a weaker monsoon would lead to lower seasonal green-up during the summer. Such changes could dampen seasonal phenology patterns, ultimately leading to the loss of a clearly detectable phenological signal over multiple years. The combination of earlier green-up and limited moisture during the initial vegetation surge would most likely lead to increasingly unfavorable conditions for native vegetation communities and propagate a conversion of native shrublands into invasive annual grasslands better able to deal with the increasingly limited moisture availability. The immediate implications of prolonged senescence and early browning further extend onto increased fire activity and the availability of easily ignitable fuel, thus raising the potential of recurrent wildfires of increasing extent and intensity across the Southwestern U.S (Westerling et al 2006, Westerling and Bryant 2007, Abatzoglou and Kolden 2011. With temperature increasingly driving contemporary drought in the Southwestern U.S (Gremer et al 2015), the interaction between precipitation and temperature is essential in determining the amount, timing and spatial extent of water available to vegetation. Observations from the Santa Rita Experimental Range have illustrated the differential effects of recent climate change, as a significant decrease in perennial grass cover and a visible reduction of vegetation density were observed due to warming temperatures and lack of essential spring and monsoon precipitation in recent years. Repeat photography has shown that the driest monsoon on record in 2020, followed by dry and hot conditions in spring 2021, led to severe vegetation die-back until the return of precipitation during the summer of 2021 resulted in a sharp increase in greenness the next year, as plants maximized their response to favorable conditions (https://cals.arizona.edu/ srer/photos.html). Similarly, our simulations have shown that protracted periods of reduced precipitation would demonstrably lead to reduced peak greenness which could go as far a as lead to early mid-summer senescence. Depending on the extent of the projected changes to precipitation timing and frequency, as well as local groundwater decline, the adaptive capabilities of dryland plant functional groups will be increasingly challenged. The usual rapid responsiveness of grasses and shrubs to favorable conditions may be compromised to a point where they are unable to respond to future precipitation events following prolonged drought periods, while at the same time riparian trees may not be able to match the speed of groundwater decline.
Based on historic climate-phenology relationships we demonstrated that the explicit inclusion of hydroclimate forcing within an empirical modeling framework is effective in capturing current and potential future vegetation responses, including seasonal changes as well as interannual variations driven by climate. A caveat to the present model is that it currently does not include other atmospheric interactions and feedbacks on phenology, such as the carbon and energy cycles. Still, the results have several implications for modeling climate-phenology interactions. The analysis provided not only insight into differential vegetation responses of various plant functional groups, but we outlined a simple empirical approach for including more dynamic phenology models into other existing ecosystem and water balance models, which might currently lack a dynamic parameterization of phenology. Especially in dryland riparian environments, improved water balance estimations have recently emerged through refined modeling frameworks (Quichimbo et al 2021), which consider all aspects of the dryland water balance, but still lacking a dynamic vegetation parameterization. The inclusion of dynamic climate-vegetation responses in such water balance models would inherently enhance the models' abilities to explore climate-soil-vegetation feedbacks on the dryland water balance as well as long-term seasonal greenness responses, and potential thresholds of vegetation mortality and succession under future climate change. Furthermore, an improved dynamic parameterization of vegetation within ecohydrological models would be essential to improve our understanding not only of historical ecohydrology but the range of impacts on climate-vegetation interactions and water balance estimates.

Conclusion
In this study we explored hydroclimate-phenology interactions for different plant functional types, which informed the development of a simple, parsimonious empirical phenology model. The proposed framework features several key features, such as a link to antecedent and concurrent climate variables, a continuous prediction of phenology throughout the season, the use of remote sensing data over an extended spatial and temporal scale, and a relative low level of complexity that can be transferred to different environments. Simple climate change scenarios have demonstrated the functionality of the model to explore vegetation responses under future climate change. While antecedent precipitation was the key driver of greenness in grasses and shrubs, depth to groundwater and to a certain extent antecedent available precipitation, were influential in controlling greenness in riparian trees. The results have several implications for climate-phenology interactions, providing not only insight into differential vegetation responses of various plant functional groups, but outlining a possible approach for including vegetation phenology models into other existing ecosystem models, which might currently lack a dynamic parameterization of phenology. To fully realize the benefits of this approach of parameterizing phenology of different plant functional groups, further model development will focus on adapting the proposed framework to be included within a dryland water partitioning model (Quichimbo et al 2021). This will allow us to explore the effects of future climate change on climate-phenology interactions and its associated impacts on the dryland water balance under varying rainfall regimes over larger spatial and temporal scales.

Data availability statement
The data that was used in this study is publicly available from different sources cited in the manuscript.
The data that support the findings of this study are openly available at the following URL/DOI: 10.5281/ zenodo.6791901.