A brown wave of riparian woodland mortality following groundwater declines during the 2012–2019 California drought

As droughts become more frequent and more severe under anthropogenic climate change, water stress due to diminished subsurface supplies may threaten the health and function of semi-arid riparian woodlands, which are assumed to be largely groundwater dependent. To better support the management of riparian woodlands under changing climatic conditions, it is essential to understand the sensitivity of riparian woodlands to depth to groundwater (DTG) across space and time. In this study, we examined six stands of riparian woodland along 28 km of the Santa Clara River in southern California. Combining remote sensing data of fractional land cover, based on spectral mixture analysis, with historical groundwater data, we assessed changes in riparian woodland health in response to DTG during the unprecedented 2012–2019 California drought. We observed a coherent ‘brown wave’ of tree mortality, characterized by decreases in healthy vegetation cover and increases in dead/woody vegetation cover, which progressed downstream through the Santa Clara River corridor between 2012 and 2016. We also found consistent, significant relationships between DTG and healthy vegetation cover, and separately between DTG and dead/woody vegetation cover, indicating that woodland health deteriorated in a predictable fashion as the water table declined at different sites and different times. Based on these findings, we conclude that the brown wave of vegetation dieback was likely caused by local changes in DTG associated with the propagation of precipitation deficits into a depleted shallow alluvial aquifer. These factors suggest that semi-arid riparian woodlands are strongly dependent on shallow groundwater availability, which is in turn sensitive to climate forcing.


Introduction
Riparian trees often rely on shallow groundwater to meet their water needs, so water table declines during increasingly frequent and severe drought conditions threaten the health and function of riparian woodlands (Diffenbaugh et al 2015, Meixner et al 2016, Rohde et al 2017, Williams et al 2020. Though many studies have examined the sensitivity of individual species to drought and local groundwater availability (e.g. Stromberg et al 1996, Singer et al 2014, Sargeant and Singer 2016, Pettit and Froend 2018, Skiadaresis et al 2021, few studies have considered the landscape-scale responses of riparian woodlands to drought across space and time, and even fewer in relation to direct groundwater measurements (but see Rohde et al 2021). These knowledge gaps are relevant to dryland watersheds around the globe, where convergent pressures on water resources from agriculture, urban development, and climate change are increasing (Taylor et al 2013, Rohde et al 2017, Rateb et al 2020.
Riparian woodlands are ecologically important plant communities that provide habitat for sensitive animal species (Kus 1998, Merritt and Bateman 2012, Bateman and Merritt 2020, promote plant biodiversity (Stromberg et al 1996, Stromberg andMerritt 2016), and contain a disproportionately large amount of the biomass in dryland watersheds (Swetnam et al 2017, Matzek et al 2018, Dybala et al 2019. They form in convergent topographic zones, which serve as hydrologic refugia and are somewhat buffered from normal climatic variability (Brooks et al 2015, Hoylman et al 2019. Riparian tree species are typically phreatophytes, which have taproot systems that extend up to 5 m down to the capillary fringe above perennial water tables (Stromberg 2013, Rohde et al 2017. Phreatophytes are extremely sensitive to water availability at all life stages (Mahoney and Rood 1998, Singer et al 2013. If their root systems lose contact with the alluvial water table, phreatophytes commonly exhibit stomatal closure, leaf abscission, branch dieback, and xylem cavitation (Scott et al 1999, Leffler et al 2000, Rood et al 2000. Prolonged water table declines, for example during extreme drought conditions, can lead to whole-plant mortality if a tree's hydraulic system cannot maintain a favorable water balance as groundwater supply declines (Scott et al 1999, Cooper et al 2003. From 2012 to 2019, California experienced the most severe drought in its paleoclimate record (Robeson 2015). Meteorological (Thomas et al 2017), and upland canopy water content (Asner et al 2016) throughout the region. While the drought is known to have generated mass die-off of upland trees (e.g. Goulden and Bales 2019), there is a notable lack of quantitative assessments of droughtinduced mortality of lowland riparian phreatophytes, particularly within dryland regions. Documenting large-scale ecological die-offs, particularly in ecosystems that are buffered from climate impacts by their hydrogeomorphic setting, is important for identifying global signals of forests being pushed past their tolerance for environmental change ( Remote sensing is a powerful tool for analyzing the sensitivity of vegetation health to changes in the water availability, but few remote sensing studies have quantified the sensitivity of groundwater-dependent ecosystems to water table declines (e.g. Barron et al 2014, Huntington et al 2016. In this study, we combined time series remote sensing imagery with data from groundwater monitoring wells to investigate the fate and trajectory of riparian woodlands in southern California during the unprecedented 2012-2019 California drought. We used spectral mixture analysis (SMA) to discriminate between healthy and dead vegetation cover in remote sensing imagery. We then analyzed the relationship between vegetation cover and depth to groundwater (DTG) in a range of woodland stands that represent a gradient of groundwater availability. These data enabled us to characterize the trajectory and the spatial progression of drought impacts across an ecologically and economically important river corridor, and to monitor the initial drought recovery in the riparian woodlands.

Study area
The Santa Clara River flows 132 km from the Mojave Desert to the Pacific Ocean and has a catchment covering 4200 km 2 in Ventura and Los Angeles Counties, California (Beller et al 2016). Mean annual precipitation ranges from 200 to 800 mm, with the wettest regions at high elevations (catchment relief 2700 m) and near the coast (Downs et al 2013). The basin has a Mediterranean climate with cool, wet winters and warm, dry summers, and many reaches of the Santa Clara River are ephemeral (i.e. flowing only part of most years). Winter rainfall produces flashy flows, and more than half of the annual discharge occurs during a small number of precipitation events (Downs et al 2013, Beller et al 2016. The river corridor has been subject to extensive urban and agricultural development over the last century, but the main stem of the river has not been severely controlled by engineering structures, making it the largest river in southern California that is mostly free flowing (

Study sites
We identified six stands of Populus-Salix riparian woodlands in the lower Santa Clara River floodplain that are thought to be supported by perennial shallow aquifers (figure 1; Beller et al 2016). The woodlands range in area from 7 to 120 ha, and they represent the most substantial woodlands that were present before the 2012-2019 drought (Beller et al 2016). The study sites are distinguished by transitions in hydrology or river management that facilitate shallow groundwater depths. Site boundaries were manually digitized in GIS software using 2012 aerial imagery acquired by the National Agricultural Imagery Program (U.S. Department of Agriculture 2012).

Depth to groundwater
We calculated DTG at each of the study sites using measurements from nearby wells acquired from the California Department of Water Resources (https://sgma.water.ca.gov/webgis/ ?appid=SGMADataViewer), United Water Conservation District, and the County of Ventura. We also used unpublished data from shallow monitoring wells that were installed by members of our team at two of the study sites (tables S1, S2; figures S1-S4 (available online at stacks.iop.org/ERL/16/084030/mmedia)). The shallow wells are ∼3 m deep and were installed between 2015 and 2020. They are manually measured twice per month. We used two different protocols to calculate DTG at the study sites, depending on the well data availability for each site (see supplementary material).

Remote sensing data acquisition and processing
SMA was used to map the fractional cover of green vegetation (GV), non-photosynthetic vegetation (NPV; i.e. dead and woody plant material), and soil in the Santa Clara River floodplain (Smith et al 1990, Roberts et al 1998. The SMA model was calibrated using data from the Airborne Visible/Infrared Imaging Spectrometer (Green et al 1998)  The SMA model generated estimates of the fractional cover of GV, NPV, and soil within each pixel for each image. While the SMA method did not classify species cover, qualitative observations of species cover were made during field visits to the sites between 2017 and 2021, and by manually examining high resolution aerial imagery captured before, during, and after the drought.

Analysis of drought effects on vegetation
The fractional cover data and the groundwater data were used to conduct two analyses. First, the GV and NPV fractions were used to examine the spatial and temporal trends of woodland mortality along the river corridor. Mortality was indicated by a decrease in GV fractions and an increase in NPV fractions (Huang et al 2019). Significant differences in land cover fractions for each study site across time were identified using a Kruskal-Wallis test and a post-hoc Dunn's test with a Holm adjustment (α = 0.05). Second, we quantified the sensitivity of GV and NPV fractions to DTG. The median GV and NPV fractions were calculated for each study site and each year. We used the DTG measurements that were closest in time to the Landsat image acquisition dates. The difference between the well measurement dates and the image acquisition dates ranged from 0 day to 36 days with a median of 13 days. The sensitivity analyses were divided into two distinct time spans. The first time span was limited to data from 2011 to 2016, representing the period when the drought became progressively more severe, as evidenced by increasing DTG, below-average soil moisture (Warter et al 2021), and decreasing Standardized Precipitation Evaporation Index (SPEI). The 2011-2016 observations were pooled across sites and years, and mixed effect logistic-binomial regression (Gelman and Hill 2007) was used to determine if DTG is a significant predictor of GV and NPV fractions. Site was included as a random effect in the models to account for local  influences on the vegetation unrelated to groundwater (see supplementary material).
The second time span was limited to data from 2017 and 2018, which represents a period of early drought recovery in the riparian woodlands. Substantial rainfall in the winter of 2016-2017 reduced DTG and increased soil moisture in the region (Warter et al 2021). As a result, the ecological drought began to subside, even though meteorological drought conditions persisted until 2019. The sensitivity of an ecological response to an environmental driver often differs based on the direction of the change. Differing sensitivities during decline and recovery phases can result in hysteresis in ecological systems (Beisner et al 2003, Andersen et al 2009. The 2017-2018 data were used to determine if there were differences in the sensitivity to DTG during the early stages of drought recovery as compared to the drought onset years of 2011-2016. The observed land cover fractions from 2017 and 2018 were compared to the values predicted by the regression model that was calibrated using data from 2011 to 2016. This indicated whether sensitivity to DTG differed during the two phases. Mean absolute error (MAE) was used to quantify the difference between the observed and predicted values during the two phases.

Groundwater level changes during and after drought
In 2011, maximum DTG at the six study sites ranged from 1.1 to 4.4 m. The DTG increased (i.e. the water table declined) at all six study sites between 2011 and 2016, but there was substantial spatial and temporal variability in DTG along the river corridor (figure 3;  table S3). The changes in DTG were mediated by the interaction between climatic forcings and basin geomorphology. The Fillmore Ciénega site sits at the boundary of the Piru and Filmore groundwater subbasins (figure 1), where the deposits of permeable alluvium become substantially narrower and shallower and force groundwater to the surface (Mann Jr 1958, Reichard et al 1999. Surface flow between the Piru and Fillmore subbasins decreased between 2011 and 2013 and stopped between 2014 and 2016 (United Water Conservation District 2017). As a result, shallow lateral recharge of the Fillmore subbasin was likely reduced or eliminated during the peak of the drought. Maximum DTG at Fillmore Ciénega and Sespe Confluence increased by 11.9 m and 12.7 m, respectively, between 2011 and 2016 as groundwater in the Fillmore subbasin was depleted. At the downstream end of the Fillmore subbasin, groundwater elevations were more stable. The East Grove site sits at the boundary of the Fillmore and Santa Paula subbasins, where constrictions in the deposits of unconsolidated alluvium again force groundwater to the surface (Reichard et al 1999). Maximum DTG at the East Grove site only increased by 0.9 m between 2011 and 2016. Surface flow at the boundary between the Fillmore and Santa Paula subbasins did not approach zero until 2016 (United Water Conservation District 2017). The Hanson, Freeman Upstream, and Freeman Downstream sites are located in the Santa Paula subbasin, where the shallow aquifer sits on top of impermeable deposits that prevent groundwater from percolating into deeper aquifers (Reichard et al 1999, Hanson et al 2003. Groundwater elevations in the Santa Paula subbasin remained relatively stable throughout the drought. At Hanson, Freeman Upstream, and Freeman Downstream, maximum DTG increased by 3.5 m, 3.0 m, and 5.0 m (respectively) between 2011 and 2016. In 2016, maximum DTG at the six sites ranged from 2.0 to 17.1 m. Substantial rainfall in the winter of 2016-2017 reversed groundwater trends and caused DTG to decrease (i.e. the water table rose) at all study sites, but some sites did not approach pre-drought DTG until 2019.

Vegetation cover change during drought
Several riparian woodlands exhibited large decreases in GV cover and large increases in NPV cover from 2011 to 2016, indicating widespread drought-induced mortality (figure 4, tables S4-S6). Fillmore Ciénega and Sespe Confluence exhibited the largest decreases in GV fractions and the largest increases in NPV fractions. Sites farther downstream were less affected by the drought and experienced smaller changes in land cover.
There was a distinct spatial pattern and temporal trend of woodland mortality that occurred both within and across study sites in the Fillmore subbasin (figure 5). Widespread mortality first occurred in 2013 at the most upstream study site, Fillmore Ciénega (figure S5). A wave of mortality then traveled 13 km west (downstream) across the Fillmore subbasin between 2013 and 2016. The wave of mortality can be seen within and across individual study sites, and it is especially distinct within the Fillmore Ciénega site. By 2015, the wave of mortality reached the area immediately upstream of the East Grove site (figure S6). By 2016, all three areas had experienced widespread mortality, and the Fillmore Ciénega and Sespe Confluence sites experienced nearcomplete mortality of their riparian woodlands. Sites in the downstream Santa Paula subbasin were relatively stable throughout the drought, and no distinct spatial pattern of woodland mortality in the Santa Paula subbasin was observed (figure S7).

Groundwater declines and plant health
DTG was significant predictor of the GV and NPV fractions for 2011-2016 (table 2). There was a significant negative relationship between DTG and the GV fractions (p < 0.001; figure 6(a)), indicating that GV decreased as the water table declined. There was also a significant positive relationship between DTG and the NPV fractions (p < 0.001; figure 6(b)), indicating that dead and woody plant cover increased as the water table declined. Taken together, these remote sensing metrics indicate leaf shedding, increased litter, exposed branches, and, in some cases, complete mortality as the drought progressed (Adams et al 1995).

Site-based differences in early drought recovery
We also compared the 2017-2018 GV fractions to the values that were predicted by the regression model, which was calibrated using data from 2011 to 2016. This revealed whether the sensitivity to DTG differed during the decline and recovery phases. The 2017-2018 GV fractions for Fillmore Ciénega deviated   S8). The other sites, where the native woodlands remained largely intact (e.g. figure S9), exhibited consistent sensitivity to DTG during the decline (MAE = 0.02-0.07) and recovery (MAE = 0.03-0.09) phases.

Discussion
During the 2012-2019 California drought, riparian woodland mortality in the Santa Clara River floodplain followed a coherent spatial pattern and temporal trend that occurred across the river corridor and mirrored the apparent trend in DTG. Mortality first occurred at the upstream side of the Fillmore Ciénega site as flow between the Piru and Fillmore subbasins decreased around 2013. A distinct 'brown wave' of mortality then traveled west (i.e. upstream to downstream) between 2013 and 2016 as groundwater in the Fillmore subbasin was progressively depleted. The brown wave stopped just upstream of the East Grove site, where groundwater elevations were relatively stable and surface flow was maintained throughout most of the drought (United Water Conservation District 2017 Groundwater serves as a crucial link in the chain of drought propagation from meteorological drying to plant responses by riparian phreatophytes. The 2012-2019 California drought was caused by record low precipitation and record high temperatures, which reduced water inputs to ecosystems and increased evaporative demand (Diffenbaugh et al 2015, Warter et al 2021. The drought generally reduced groundwater recharge (Harlow and Hagedorn 2018) and caused groundwater elevations to decline, but subsurface water fluxes are spatially variable and are mediated by several factors including agricultural water use and runoff, river regulation, soil texture, and bedrock geology. During drought conditions, the hydrological drivers of groundwater elevation interact with meteorological trends to produce distinct spatial patterns and temporal trends of groundwater change (Jencso andMcGlynn 2011, Harlow andHagedorn 2018).
The spatiotemporal variability in groundwater elevation caused varying physiological responses in the riparian woodlands in the Santa Clara River floodplain. When water tables decline by several meters, the root systems of riparian trees can lose access to groundwater (Stromberg 2013). Riparian tree species are poorly adapted to drought and can experience catastrophic xylem cavitation at relatively high (i.e. close to zero) vapor pressure deficits (Fichot et al 2015). To mitigate and prevent this often-irreversible change, trees undergo a series of physiological changes to maintain a favorable water balance in the face of declining water supply (Rood et al 2003). Within minutes, they can regulate stomatal conductance to limit transpirational water loss (Horton et al 2001, Amlin and Rood 2003, Pivovaroff et al 2018. Leaf abscission and branch dieback, though detrimental to woody plants in the short term, can serve as long-term survival strategies that help trees reduce water demand and prevent the loss of xylem water conductance (Scott et al 1999, Rood et al 2000, Cooper et al 2003. These physiological responses may not be adequate to mitigate water stress from large and sudden water table declines, which often cause hydraulic failure and whole-plant mortality (Scott et al 1999, Lite and Stromberg 2005, Tai et al 2018. The brown wave is likely an emergent property of individual plants responding to localized changes in DTG, as evidenced by the fine-scale changes in plant health and the strong statistical relationships between the land cover fractions and DTG. Few studies have examined the spatial evolution of riparian woodland responses to water table declines (but see Stromberg et al 1996, Scott et al 1999, Tai et al 2018. The current understanding of phreatophyte sensitivity to DTG is largely derived from field measurements (e.g. Horton et al 2001, Rood et al 2011), laboratory experiments (e.g. Leffler et al 2000, , and models based on phreatophyte physiology (e.g. Tai et al 2018). These data have resulted in a conceptual model that suggests that there is a highly non-linear relationship between DTG and plant health, whereby plant health degrades rapidly when DTG increases beyond some critical threshold (e.g. Shafroth et al 2000, Horton et al 2001, Lite and Stromberg 2005. In contrast with previous studies, our observations indicate that there is a mostly linear relationship between DTG and plant health at a stand scale when DTG is less than 10 m (i.e. figure 6(c)). The difference in the shape of the observed relationships may indicate a scale dependence of the analysis. While many field and laboratory-based studies examine individual plants belonging to selected species, remote sensing data detects many plants belonging to many species within a pixel (Kibler et al 2019). The observed plants that compose riparian woodlands likely have varying structures, life histories, and tolerances for groundwater decline (Stromberg et al 1996, Stromberg andMerritt 2016). Nonetheless, the observed sensitivity to absolute DTG and DTG change (table S3) were generally consistent with previously reported values. At Fillmore Ciénega and Sespe Confluence, vegetation cover stopped changing as a function of DTG between 2015 and 2016, which may indicate a fundamental limit beyond which phreatophytes experience complete mortality and the health of the surviving non-phreatophytic vegetation becomes decoupled from DTG. The recovery of phreatophytes at these sites depends on the water table returning to shallow depths and seed sources for the germination of new seedlings (Stella et al 2006).

Conclusion
The statistical analyses presented here provide some of the first robust estimates of the sensitivity of riparian woodlands to DTG. Our findings also reveal that DTG trends can be highly variable during extreme drought conditions, even within the same river corridor, which can result in distinct spatial patterns and temporal trends of plant mortality in riparian woodlands. Quantifying the sensitivity of riparian ecosystems to groundwater change will become increasingly important as anthropogenic climate change increases the frequency and severity of drought conditions across the western United States (Diffenbaugh et al 2015, Rohde et al 2017, Williams et al 2020. The widespread mortality observed during the brown wave mirrors the dynamics of mass die-offs that have occurred in upland forest ecosystems (e.g. Allen et al 2010, Goulden and Bales 2019). Anthropogenic climate change, shifts in water availability, and other environmental forcings are overwhelming the resiliency of ecosystems that are typically buffered from climatic variability (Allen et al 2015). Quantifying the sensitivity of both upland and lowland forests to hydroclimatic change will improve our ability to predict critical shifts in ecosystem structure and function in the coming decades.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.