Drought resistance across California ecosystems: evaluating changes in carbon dynamics using satellite imagery

Drought is a global issue that is exacerbated by climate change and increasing anthropogenic water demands. The recent occurrence of drought in California provides an important opportunity to examine drought response across ecosystem classes (forests, shrublands, grasslands, and wetlands), which is essential to understand how climate influences ecosystem structure and function. We quantified ecosystem resistance to drought by comparing changes in satellitederived estimates of wateruse efficiency (WUE = net primary productivity [NPP]/evapotranspiration [ET]) under normal (i.e., baseline) and drought conditions (ΔWUE = WUE2014 − baseline WUE). With this method, areas with increasing WUE under drought conditions are considered more resilient than systems with declining WUE. Baseline WUE varied across California (0.08 to 3.85 g C/mm H2O) and WUE generally increased under severe drought conditions in 2014. Strong correlations between ΔWUE, precipitation, and leaf area index (LAI) indicate that ecosystems with a lower average LAI (i.e., grasslands) also had greater Cuptake rates when water was limiting and higher rates of carbonuptake efficiency (CUE = NPP/LAI) under drought conditions. We also found that systems with a baseline WUE ≤ 0.4 exhibited a decline in WUE under drought conditions, suggesting that a baseline WUE ≤ 0.4 might be indicative of low drought resistance. Drought severity, precipitation, and WUE were identified as important drivers of shifts in ecosystem classes over the study period. These findings have important implications for understanding climate change effects on primary productivity and C sequestration across ecosystems and how this may influence ecosystem resistance in


INTRODUCTION
Drought affects ecological systems across every climatic zone worldwide and is exacerbated by climate change and increasing anthropogenic water demands (Mishra and Singh 2010). Characterized by below-normal precipitation (Dai 2011), meteorological drought results from complex interactions between the atmosphere and hydrologic processes within the biosphere. Unlike aridity, which is a permanent feature of climate (Wilhite 1992), drought is a temporary extreme event (Palmer 1965, Mishra andSingh 2010) that can persist for extended time periods (months to years; Mishra and Singh 2010). Drought can cause significant changes in ecosystem productivity and water dynamics, and it is one of the most economically and ecologically disruptive extreme events affecting millions of people globally (Dai 2011).
In California, the most recent drought began in 2012, and during the summer of 2014, ~80% of the state was in extreme to extraordinary drought and ~100% was in severe drought or worse (U.S. Drought Monitor). Combined with the diversity of natural ecosystems, multiple years of extended severe drought and the recent occurrence of the most extreme droughts on record (Diffenbaugh et al. 2015) make California an important case study to examine variations in drought resistance across ecosystems.
Here, ecosystem resistance is the capacity to absorb disturbance (i.e., drought) and retain the same function (i.e., productivity) and sensitivity (i.e., water-use efficiency [WUE]; Angeler and Allen 2016). WUE links the biological (i.e., photosynthesis and transpiration) and physical (i.e., evaporation) processes that control carbon and water dynamics, and is defined here as net primary productivity (NPP; g C/m 2 ) per amount of water lost (evapotranspiration: ET; mm/m 2 ). Drought suppresses both carbon and water dynamics simultaneously (Ponce-Campos et al. 2013, Yang et al. 2016. However, the sensitivity of the different biological and/or physical processes that influence productivity and ET depends on ecosystem type and other confounding environmental factors (Lu and Zhuang 2010, Zhu et al. 2011, Tang et al. 2014, Yang et al. 2016. Across ecosystems, WUE generally changes with precipitation (Huxman et al. 2004, Emmerich 2007, Tian et al. 2010, Yang et al. 2016. Often equated in a simplistic manner with drought resistance (Blum 2005), a high WUE translates to a greater capacity to maintain productivity under stress (Blum 2005, Ponce-Campos et al. 2013. In response to changes in conditions, WUE increases with aridity (Huxman et al. 2004, Reichstein et al. 2007, Bai et al. 2008, Lu and Zhuang 2010, Zhu et al. 2011, Ponce-Campos et al. 2013, Yang et al. 2016, and if drought becomes severe enough, a breakdown in ecosystem resistance can lead to a reduction in WUE (Reichstein et al. 2002, Lu and Zhuang 2010, Zhu et al. 2011, Yang et al. 2016) and ecosystem type conversions (Yang et al. 2016). An interesting measure of ecosystem functionality (Emmerich 2007), evaluating shifts in WUE over time under nondrought and drought conditions can provide a good approximation of ecosystem resistance to drought.
We evaluate drought resistance across California ecosystem classes (forest, shrubland, grassland, and wetland ecosystems) over 12 years (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) by quantifying deviations in WUE in 2014 from WUE under normal climate conditions. Spatial dynamics and interannual variability in WUE at large scales have rarely been quantified (Lu and Zhuang 2010, Zhu et al. 2011, Ponce-Campos et al. 2013, Tang et al. 2014, Huang et al. 2015 due to complex interactions between water and C and uncertainty in the interactive effects of multiple environmental factors on WUE (Tian et al. 2010). Here, we evaluate changes in satellite-derived WUE in response to drought to measure drought resistance across California ecosystems. We hypothesize that drought resistance in California will have a positive relationship with WUE under normal climate conditions; namely, that ecosystems with a high WUE under normal climate conditions will be most resistant to drought. Studies have shown that WUE increases with drought severity (Ponce-Campos et al. 2013) and that waterlimited ecosystems have higher WUE (Huxman et al. 2004, Reichstein et al. 2007, Ponce-Campos et al. 2013. We aim to (1) quantify spatiotemporal patterns in drought severity using the self-calibrating Palmer Drought Severity Index (scPDSI), (2) monitor drought resistance using changes in WUE, and (3)  Climate change projections indicate that extreme events will become more common in the future (IPCC 2013), making it important that we understand how ecosystems respond to these events and the potential feedbacks to radiative forcing. A critical link between C and water cycles in terrestrial ecosystems, WUE has been identified as an effective way of assessing ecosystem response to climate change (Baldocchi 1994, Hu et al. 2008, Kuglitsch et al. 2008, Beer et al. 2009, Niu et al. 2011. To our knowledge, this is the first study to examine shifts in WUE using satellite imagery and relate them to drought across ecoregions (Bailey 1995) and ecosystem classes in the state of California. This research is essential to enhance our understanding of ecosystem response to drought and how carbon dynamics change with major shifts in climate and extreme events. An analysis of severe drought effects on ecosystem function is almost completely lacking, limiting our understanding of drought resistance and how it changes with increasing drought severity. Evaluating the hydroclimatic thresholds that reduce ecosystem resistance will improve our ability to predict the consequences of increasing aridity. The findings of this study are relevant to California and more broadly to other mediterranean ecosystems around the world, which face increasing threats from drought.

Study site
California is home to a diversity of ecosystems ( Fig. 1), where ecosystem structure, function, and C dynamics are driven by differences in hydroclimate, topography, and land use. Within the humid temperate domain, California stretches across 20 ecoregions that span the mediterranean division (Bailey 1995). This mediterranean division is subject to wet winters and dry summers that often contain 2-4 months without precipitation. Drought is a natural occurrence in California, and ecosystems are likely to exhibit varying levels of drought resistance. Shrublands (46%) account for the greatest portion of natural area in California followed by forests (43%), grasslands (11%), and finally wetlands (<1%; Fig. 1). From 2002 to 2013, shifts in ecosystem classes occurred over ~ 32% of the study area.

Ecosystem class
We used the Moderate Resolution Spectroradiometer (MODIS) MCD12Q1 land cover type data to identify forests, shrublands, grasslands, and wetlands in California (Appendix S1: Tables S1 and S2; Fig. 1a). The most recent annual land cover data available (2012) defined the ecosystem class for the study, and we used annual land cover data (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) to evaluate changes in ecosystem class. The MODIS land cover type product is produced using an ensemble-supervised classification algorithm with techniques to stabilize classification results across years to reduce variation not associated with land cover change (Friedl et al. 2010). The classification algorithm includes spectral and temporal information from MODIS bands 1-7 (Huete et al. 2002) supplemented by the enhanced vegetation index and MODIS land surface temperature (Friedl et al. 2010). Year-to-year variability in phenology and disturbances such as fire, drought, and insect infestations leads to high variability that is difficult to consistently characterize the spectral signature of ecosystem classes. These effects make it harder to discern classes that are ecologically proximate and arise from poor spectral-temporal separability in MODIS data (e.g., mixed forest and deciduous broadleaf forest). To address this, the MCD12Q algorithm imposes constraints on year-to-year variation in classification results at each pixel using posterior probabilities associated with the primary label in each year (Friedl et al. 2010). If the classifier predicts a different class from the preceding year, the class label is changed only if the posterior probability associated with the new label is higher than the probability associated with the previous label (Friedl et al. 2010). To avoid propagating incorrect or out-of-date labels in areas of change, a threeyear window is used.
Classification errors are largely concentrated among classes that encompass ecological and biophysical gradients (Friedl et al. 2010). In this study, we aggregated classes into major ecosystem types by reclassifying natural ecosystems into four classes (i.e., forests, shrublands, grasslands, and wetlands; Appendix S1: Table S2). We excluded areas that were beyond the scope of this study (i.e., urban, croplands, waterbodies, and snow). Because water subsidies in agricultural systems would distort drought effects, we also excluded all areas classified as crops from this analysis. A detailed description of the MODIS MCD12Q1 product can be found at https://lpdaac.usgs.gov/ dataset_discovery/modis/modis_products_table/ mcd12q1 (Friedl et al. 2010).

Drought indices and climate conditions
We used remotely sensed drought and climate indices to quantify drought effects across the California landscape. To evaluate changes in drought condition across ecosystem classes and over time, we used monthly scPDSI data, made available through the Western Regional Climate Center (Appendix S1: Table S1). The scPDSI ranges from −5.0 to 5.0. The number shows the magnitude and the sign denotes (+) wetter than average or (−) drier than average conditions for a location based on historical climate and sensitivity to changes in water availability (Wells et al. 2004). Values of scPDSI between −0.4 and 0.4 denote average conditions and absolute values greater than 4 represent extreme conditions. Unlike earlier versions of the PDSI, extreme conditions occur based on the history of the location and are not determined relative to a default location (Wells et al. 2004). The scPDSI allows more exact comparisons between location and times and is a more accurate index compared with PDSI for extreme events (Wells et al. 2004). We also obtained PRISM precipitation (annual total mm) data sets from the PRISM Climate Group (Appendix S1: Table S1). Shrublands (46%) account for the greatest portion of natural area in California followed by forests (43%), grasslands (11%), and finally wetlands (<1%). White areas represent land cover classes that were not included as part of this study (i.e., urban, croplands, waterbodies, snow, and classes that changed 2002-2012 A detailed description of the PRISM precipitation product can be found at http://prism.nacse.org/ documents/PRISM_datasets_aug2013.pdf.

Ecosystem measures of C, ET, and LAI
To estimate ecosystem-level measures of C dynamics, we accessed MODIS MOD17 data that were processed and made available by the Numerical Terradynamic Simulation Group at the University of Montana (Appendix S1: Table  S1). Terrestrial NPP quantifies the amount of atmospheric C fixed by plants and accumulated as biomass. MOD17 estimates of NPP are based on the radiation-use efficiency logic suggested by Monteith (1972). The algorithm uses a radiation conversion efficiency concept to estimate gross primary productivity (GPP) from satellitederived FPAR (from MOD15: NDVI = FPAR) and independent estimates of photosynthetically active radiation (PAR) and other surface meteorological fields from NASA Global Modeling Assimilation Office climate data (temperature and vapor pressure deficit; Eq. 1), and the subsequent estimation of growth (R g ) and maintenance (R m ) respiration terms that are subtracted from GPP to arrive at annual NPP. ε max is the PAR conversion efficiency and translates FPAR, the fraction of PAR absorbed by the surface, to biomass produced. Radiation-use efficiency parameters (ε max and T min and VPD scalars) are extracted from the Biome Properties Look-Up Table (BPLUT) and are based on the land cover product (MOD12; Appendix S1: Table  S2). VPD is the only variable in the algorithm directly related to environmental water stress, and studies have shown that the VPD-based water stress estimate in MOD17 is adequate to explain the magnitude and variability of water stress (Mu et al. 2007).
Daily net photosynthesis (PSN net ) subtracts leaf and fine root respiration from GPP (Eq. 2).
NPP is the annual sum of daily PSN net minus the cost of growth (R g ) and maintenance (R m ) of living cells in permanent woody tissue (Eq. 3).
The maintenance respiration (MR) and growth respiration (GR) components arise from allometric relationships linking daily biomass and annual growth of plant tissues to satellite-derived estimates of leaf area index (LAI, MOD15). Satellitebased observations of NPP provide a quantitative measure of spatial patterns and seasonal to interannual variability in vegetation activity (Heinsch et al. 2006). MOD17 products have been validated at the site level using a number of eddy covariance flux tower measurements across different climatic regimes and biome types (Turner et al. 2003, 2006, Heinsch et al. 2006. A detailed overview of the MOD17 algorithm, quality control, and filling missing data can be found in the MOD17 User's Guide (Running et al. 2004, Zhao et al. 2005, Running and Zhao 2015. As an estimate of ET (mm m −2 yr −1 ), we used the MODIS global ET product (MOD16) that is also a part of the NASA Earth Observing System (EOS) project and made available by the Numerical Terradynamic Simulation Group at the University of Montana (Appendix S1: Table  S1). A vital component of the water cycle (Mu et al. 2007), ET is the sum of transpiration (linked to GPP) and evaporation from wet vegetation and soil surfaces (not linked to GPP; Kuglitsch et al. 2008). Spatial variations in precipitation and ET are a critical component of drought detection and assessment (McVicar andJupp 1998, Mu et al. 2007). Details on the processing of MOD16, quality control, and filling missing data can be found in Mu et al. (2007).
Estimates of MOD15 LAI were obtained from the NASA/EOS project (Appendix S1: Table S1). Derived via radiative transfer methods (Myneni et al. 2003, Heinsch et al. 2006, LAI is a dimensionless ratio of leaf area covering a unit of ground area (m 2 /m 2 ) and is the biomass equivalent of FPAR. The MOD15 LAI product is an essential input for the MOD16 (ET) and MOD17 (GPP/NPP) algorithms. Details on the processing of MOD15, quality control, and filling missing data can be found in Mu et al. (2007).
We compare changes in NPP, ET, LAI, and climatic conditions to understand drought resistance. NPP is sensitive to a number of controls (i.e., climate, plant characteristics, and disturbance) that influence each other and that are highly correlated (Field et al. 1995). As a result, MODIS products contain similar parameters in their algorithms (LAI/FPAR) that might make them highly correlated. Because they have been validated independently to show agreement with the observed ground measurements (Appendix S1: Table S1), MODIS products can still be used to assess ecosystem drought response to show how changes in LAI and ET relative to NPP are driving drought resistance. MODIS (land cover, NPP, and ET) and PRISM (precipitation) data sets were obtained from 2002 to 2014. Although all other data sets are available since 2000, MOD15 LAI data are only available starting in 2002, and thus, we picked 2002 as our starting year. We resampled all spatial data to match the resolution of the scPDSI data set (4 km) using bilinear methods of interpolation and re-projected data sets to UTM WGS84 using the R library raster (Hijmans 2015).

Water-use and carbon-uptake efficiency
Plant gas exchange is a key process shaping global hydrologic and C cycles and is often characterized by plant WUE (Assouline andOr 2013, Keenan et al. 2013). Allowing us to relate productivity to water dynamics (Webb et al. 1978), we estimate WUE annually using NPP (g C m −2 yr −1 ) and ET (mm −1 yr −1 or kg −1 yr −1 ; Eq. 1). Ecologists commonly use the ratio of the main ecosystem fluxes such as NPP (Roupsard et al. 2009), net ecosystem productivity/exchange (NEP/NEE; Monson et al. 2010, Niu et al. 2011), or GPP (Kuglitsch et al. 2008 to water loss (ET or transpiration; Law et al. 2002, Reichstein et al. 2002, Kuglitsch et al. 2008) as a measure of WUE (Eq. 4).
WUE has been recognized as an important characteristic of productivity in various natural scientific disciplines and has been used recently at the ecosystem level (Hu et al. 2008, Kuglitsch et al. 2008, Monson et al. 2010, Niu et al. 2011, Tang et al. 2014. Originally, WUE was used in leaf-and plant-scale studies and the theory, developed at the leaf scale, formed the foundation of most ecosystem and global scale models of WUE (e.g., Bonan 1996, Sellers et al. 1996, Pyles et al. 2000, Monson et al. 2010. MODIS estimates of WUE have been validated with tower-based estimates across ecosystem classes showing that observed WUE is consistent with MODIS-derived estimates of WUE (R 2 = 0.74-0.96; Tang et al. 2014).
We also evaluated fluctuations in ecosystem C-uptake efficiency (CUE), which is estimated using NPP (g C m −2 yr −1 ) and LAI (Eq. 5) and used only to show how NPP is changing with LAI.
Although the influence of LAI is incorporated into the estimate of productivity via FPAR, looking at changes in NPP and LAI estimates together we can further evaluate patterns in ecosystem response to drought by examining how NPP varies with LAI over time in addition to changes in WUE. Drought can cause declines in LAI, as plants can shed biomass in response to declining resources (Chapin et al. 2002). Thus, CUE expressed as a function of LAI allows us to monitor relative changes in ecosystem productivity with shifts in water availability.

MOD16 (ET) and MOD17 validation (NPP and GPP)
To understand how MOD16 (ET) and MOD17 (NPP and GPP) data compare with measured GPP and ET in California, we used eddy covariance data collected at 10 tower sites across California (Goulden et al. 2012). Tower sites were located in the central Sierra Nevada and San Jacinto Mountains and extend across forest, grassland, and shrubland ecosystems (Table 1). Eddy covariance data were collected, processed, and made available by the Goulden Lab at the University of California Irvine (http://faculty.sites.uci.edu/ mgoulden). Data processing information can be found in Goulden et al. (2006Goulden et al. ( , 2012. Linear comparisons between MODIS products and towerbased measures of GPP and ET were documented to show how well MODIS products represent measured conditions. We also compared MODISderived WUE to tower WUE. For the purpose of this analysis only, WUE is calculated with GPP. Throughout the remainder of the study, we calculate WUE with NPP, which includes the effects of growth and maintenance respiration.

Study design
We defined mean annual baseline values for precipitation, NPP, ET, LAI, WUE, and CUE under normal conditions (i.e., scPDSI values  We used spatial regression model estimation to determine drivers of ΔWUE using the cells with baseline data. Using a general-to-specific approach, we began with a simple linear model (using the lm function in base R) that was extended depending on the results of the Moran's I test for residual spatial autocorrelation, with the lm.morantest function in the R package spdep (Bivand and Piras 2015). If spatial autocorrelation was present, the Lagrange multiplier test statistic for spatial autocorrelation was used, applying the lm.LMtest function in spdep, to determine whether a spatial error model or a spatial lag model (SLM) was appropriate. Maximum-likelihood (ML) estimation of SLMs was carried out with the lagsarlm function in the R package spdep. We used the significance of ρ, the likelihood ratio test, and measures of fit (AIC) to evaluate SLMs. Maximum-likelihood estimation of the spatial error model (SEM) is similar to the lag procedure and implemented in the errorsarlm function in the R package spdep. We used the significance of λ, the likelihood ratio test, and measures of fit (AIC) to evaluate SEMs. An important specification test in the SEM is a test on the spatial common factor hypothesis. This exploits the property that a SEM can also be specified in spatial lag form, with the spatially lagged explanatory variables included, but with constraints on the parameters. The spatial lag form of the error model is also referred to as the spatial Durbin specification. A test on the common factor hypothesis was performed using the LR.sarlm function. This is a simple likelihood ratio test, comparing two objects for which a logLINK function exists. Explanatory variables (CUE, precipitation, NPP, LAI, and scPDSI) were strategically selected to reduce multicollinearity, and model selection was based on measures of fit (AIC). To remove the effects of possible changes in vegetation classes or errors in classification, cells that did not maintain the same ecosystem class over the study period were removed from this analysis of drought resistance and were evaluated separately to explore drivers of ecosystem class type conversions. Data used in this analysis has been archived on the Knowledge Network for Biocomplexity (Malone 2015).
Covariance between type conversions from one ecosystem class to another and drought resistance was explored using logistic regression techniques with a data set that included all areas that changed ecosystem classes over the study period. We utilized logistic regression methods to determine the probability of converting to a forest, shrubland, or grassland ecosystem on a pixel level using the function glm in base R (R Core Team 2014). Each pixel was treated as an observation, responses were coded as 1 or 0 for "success" or "failure," and maximum-likelihood methods were used for parameter estimation. With this approach, parameters were estimated iteratively until parameters that maximized the log of the likelihood were obtained. To determine whether the model displayed lack of fit, we used the ratio of the Pearson chi-square statistic to its degrees of freedom. Values closer to 1 indicated that the models fit the data well. Because raster data are comprised of adjacent pixels, the assumption of independence among observations was likely violated due to spatial autocorrelation. We incorporated latitude and longitude in models to account for both spatial autocorrelation and to examine spatial trends.
All parameters (longitude, latitude, scPDSI, precipitation, Δprecipitation, and ΔWUE) were kept in the final model for comparison.

Baseline conditions
Ecosystem classes extend across a range of conditions in California leading to complex relationships between baseline layers for ecosystem classes and ecoregions. Here, we show both the general trends in covariance between layers and compare them across and between ecosystems classes. Although WUE was positively correlated with CUE, NPP, and LAI, average baseline WUE was highest in grassland systems, which also had lower rates of mean annual NPP, precipitation, and LAI. Under baseline conditions, WUE was negatively correlated with precipitation across all ecosystem classes and ecoregions. Although mean annual baseline WUE was highest for grasslands followed by forests and shrublands, precipitation was highest for forest and lowest  Fig. 2). Forests were more likely to occur where precipitation exceeded 839 mm/yr, shrublands within the range of 360-839 mm/yr, and grasslands dominated where precipitation was less than 359 mm/yr (Fig. 2). Compared with shrublands and grasslands, forest had the highest average annual precipitation, NPP, and LAI (Table 2). CUE exhibited patterns similar to WUE and was highest in ecosystems with the lowest annual rates of precipitation and LAI (Fig. 3). Although CUE was positively correlated with NPP, grasslands had the greatest mean annual CUE followed by shrublands and forests.

Temporal and spatial patterns in drought
Based on the mean annual scPDSI for the entire state, from 1896 to 2014 (119 yr), 51 yr (43%) was in drought (scPDSI < −1). While drought is a regular occurrence in California ecosystems (Fig. 4), extreme drought represents just 2% of the sampled period (1896 to 2014). Spatial patterns in scPDSI show that drought conditions increased in severity with declining latitude (southward) and were greatest in the central portion of California (Fig. 4c). In 2014, the lowest values for scPDSI were between 34° N and 38° N, which corresponded to areas dominated by shrublands and grasslands (Figs. 1 and 4c). Capturing relative trends between ecoregions, scPDSI showed that the southern portion of the state was drier (deserts) than the cool moist forests to the north (Figs. 1 and 4c).

Quantifying drought resistance with ΔWUE
To measure drought resistance, we compared baseline WUE to WUE under extreme drought conditions in 2014 (i.e., ΔWUE; Figs. 5a and 6).

Changes in ecosystem class
From 2002 to 2012, forests (3%), shrublands (17%), grasslands (35%), and wetlands (22%) all experienced changes in ecosystem classes. Forests were most likely to convert to shrublands (2%), while shrublands, grasslands, and wetlands were more likely to convert to forests (12, 20, and 22%, respectively; Table 4). The probability of converting to a forest was positively related to precipitation and longitude; deviations in precipitation and ΔWUE decreased the probability of converting to a forest (Table 5). Precipitation and ΔWUE increased the probability of converting to a shrubland, while longitude, scPDSI, and Δpreci pitation reduced the likelihood ( Table 5). The probability of converting to a grassland increased with longitude, scPDSI, and Δprecipitation while declining with latitude, precipitation, and ΔWUE (Table 5). It is important to note that errors in classification, although reduced by aggregating MODIS original land cover classes, cannot be ruled out as influencing changes in ecosystem class.

ΔWUE as an indicator of drought resistance
Water limits plant growth and productivity in nearly all terrestrial ecosystems, and biomes differ in their sensitivity to variations in water availability (Huxman et al. 2004). Major ecosystem classes (i.e., forests, shrublands, grasslands, and wetlands) with typical dominant plant functional types are often characterized by diverging water availability (Woodward 1987 , Fig. 3). Physiological variations in leaf gas exchange and differences in environmental conditions influence ecosystem WUE (Farquhar et al. 1989, Kuglitsch et al. 2008. Our results support these findings and suggest that in California baseline WUE varies across ecosystems (0.08-3.85 g C/kg H 2 O) and was greatest in grasslands followed by forests and shrublands. WUE in California ecosystems fall within the range of values reported in other studies (Table 6).  Distinct differences in WUE and ΔWUE among ecosystems might be determined by both inherent biological characteristics and external environmental conditions.
We hypothesized that drought resistance in California would decline with decreasing baseline WUE. In general, WUE is expected to increase with increasing aridity (Huxman et al. 2004, Reichstein et al. 2007, Bai et al. 2008, Ponce-Campos et al. 2013) and a high WUE often means greater capacity to maintain productivity under stress (Blum 2005, Ponce-Campos et al. 2013. Our results indicate that for the majority of the study area and across all ecosystem classes, WUE and ΔWUE increased with aridity. A higher ET sensitivity to changes in hydroclimatic conditions is the primary reason for the negative correlation between WUE and drought severity in drought-resistant ecosystems (Yang et al. 2016). Where baseline WUE was ≤0.4, WUE declined under extreme drought conditions, suggesting that a threshold for baseline WUE exists where drought resistance is low. A baseline WUE above 2.7 is indicative of systems with very high drought resistance in California. Our results suggest that higher resistance (i.e., greater ΔWUE) is associated with lower LAI and greater C-uptake rates when water is limiting by either increasing WUE or increasing CUE. Across all ecosystem classes, WUE hints at adaptation to the local environment (Kuglitsch et al. 2008).

Ecosystem resistance to drought
In California, ecosystem resistance is lower in high productivity areas that are not adapted to maintain or increase CUE as water resources become limiting. Studies have shown that ecosystem productivity generally increases with precipitation (Paruelo et al. 1999, Huxman et al. 2004, Bai et al. 2008, Tian et al. 2010) and NPP and WUE closely match temporal and spatial precipitation patterns in NPP (Tian et al. 2010, Ponce-Campos et al. 2013. In our study, WUE was highest in the southern portion of California where drought conditions were more severe (Fig. 5). Comparing baseline WUE with WUE during extreme drought supports previous findings of an increasing trend in WUE with water stress (Ponce-Campos et al. 2013). Drier regions had higher WUE and ΔWUE (i.e., higher drought resistance) and lower LAI, GPP, and NPP. Similar to the findings of other studies (Lauenroth et al. 2000, Reichstein et al. 2007, Li et al. 2008, Niu et al. 2011, Ponce-Campos et al. 2013, our results suggest that WUE increased with decreasing precipitation, a response primarily controlled by C rather than water processes (Niu et al. 2011). Consistent correspondence between WUE and NPP (Hu et al. 2008, Niu et al. 2011 suggests that photosynthetic processes are the dominant regulator of seasonal variations in WUE, which may be responding to the effect of LAI on C assimilation (Hu et al. 2008) and differences in CUE. LAI has been shown to be significantly correlated with WUE, driving seasonal, interannual, and site-to-site variation in WUE (Hu et al. 2008) although its effect varies with GPP (temporally and spatially; Hu et al. 2008). The positive relationship observed between scPDSI and LAI suggests that reducing or maintaining a low leaf area is a major mechanism for moderating water use, especially under drought stress (Mitchell et al. 1998, Blum 2005, and ecosystems with low LAI may have adapted to drought by altering CUE as water becomes limiting. Ecosystem LAI is determined largely by the availability of soil resources (mainly water and nutrients), climate, and time since disturbance and is both a cause and a consequence of differences in NPP (Chapin et al. 2002). Drought causes plants to reduce leaf area below the potential LAI (Chapin et al. 2002), as plants often shed biomass whenever resources decline below some threshold needed for maintenance (Chapin et al. 2002). Across all ecosystem classes, our results indicate that CUE is positively correlated with WUE. Although WUE increased with aridity for the majority of the study area, declines in WUE under drought conditions were also observed. Declining WUE is indicative of a breakdown in ecosystem resistance to drought (Lu and Zhuang 2010). As drought resistance declines, GPP dominates WUE variability (Yang et al. 2016). Studies have shown that under extreme drought conditions, WUE drops below the values observed under "normal conditions" (Lu and Zhuang 2010). This change in drought resistance is most likely due to reductions in carbon assimilation capacities (i.e., impairment of electron transport and carboxylation capacity) and changes in active leaf area (Harley and Tenhunen 1991, Reichstein et al. 2002. We observed negative ΔWUE under a range of drought conditions, suggesting that the length of drought might also influence this decline in drought resistance.

Implications for climate change
As a major feedback to global climate, terrestrial carbon cycling directly affects the concentration of atmospheric CO 2 (Chapin et al. 2002). Positive feedbacks to drought occur when decreased precipitation leads to moisture depletion within ecosystems. Early signs of climate change have been observed in California where trends toward warmer temperatures and shifts in precipitation patterns have been observed (Cayan et al. 2001, Mote et al. 2005, Knowles et al. 2006, Stewart et al. 2010. Climate change simulations suggest that possible future climate scenarios include higher temperatures (+1.5°C to +4.5°C) in California by 2100 (Cayan et al. 2008) and show that warming will be greater in the summer compared with winter months with relatively small (less than ~10%) changes in overall precipitation (Cayan et al. 2008). These projections suggest that water demands will increase into the future (Diffenbaugh et al. 2015), with little to no change in precipitation. The impact of future climate on ecosystem C and water dynamics is uncertain, making it essential to evaluate drought resistance across ecosystems to understand how climate change will alter ecosystem structure and function. Our research suggests that as water becomes more limiting, ecosystem function across the landscape is at risk. Shifts in ecosystem classes have already occurred over ~32% of California, and our data suggest that this trend is influenced by precipitation and drought conditions. Changes in the distribution of wetlands across the California landscape may have already occurred.

Research limitations
This research evaluates drought effects across California ecosystems using remote sensing imagery to address hypotheses involving characteristics that enhance drought resistance. Aside from limitations imposed by using satellite imagery (e.g., aggregation at the pixel scale), inherent limitations of this study include issues addressing spatial patterns without expansive groundvalidated data sets. This research is further restricted by the lack of information on the relationship between water demands and C uptake. Without a clear understanding of water demands, we do not actually know at what point systems are water limited. Several caveats should be noted. First, the assumption is that error in the use of the MODIS parameter as a proxy for productivity is invariant across ecosystem types. Second, we are also limited by our inability to distinguish C3-and C4-dominated areas, which we would expect to show differing sensitivities to water stress. Despite these limitations, this is the first study to show drought resistance across ecosystem types in California using remote sensing-derived WUE and to show that previously observed patterns in WUE are expressed at the ecosystem level.

CONCLUSIONS
We examine changes in ecosystem WUE in response to drought to enhance our understanding of ecosystem resistance under extreme drought conditions. This study demonstrates how remote sensing products and open data can be used to monitor ecosystem resistance to drought. Our analysis establishes a framework to utilize these products to understand differences in ecosystem response to drought over large spatial scales. Furthermore, ΔWUE functions as an indicator of ecosystem resistance in light of projected climate change. ΔWUE is an important indicator of ecosystem vulnerability by signaling severe changes in resource availability that cause shifts in the structure and function of ecosystems across the landscape. Our results, along with findings from other studies, show how precipitation, LAI, CUE, and WUE play an important role in regulating C uptake and sequestration. In California, arid ecosystems had greater drought resistance compared with high productivity ecosystems in the northern portion of the state. These findings have important implications for understanding climate change effects on primary productivity and C sequestration across ecosystems and how this may influence ecosystem resistance in the future.