Quantifying Drought Sensitivity of Mediterranean Climate Vegetation to Recent Warming: A Case Study in Southern California

: A combination of drought and high temperatures (“global-change-type drought”) is projected to become increasingly common in Mediterranean climate regions. Recently, Southern California has experienced record-breaking high temperatures coupled with significant precipitation deficits, which provides opportunities to investigate the impacts of high temperatures on the drought sensitivity of Mediterranean climate vegetation. Responses of different vegetation types to drought are quantified using the Moderate Resolution Imaging Spectroradiometer (MODIS) data for the period 2000–2017. The contrasting responses of the vegetation types to drought are captured by the correlation and regression coefficients between Normalized Difference Vegetation Index (NDVI) anomalies and the Palmer Drought Severity Index (PDSI). A novel bootstrapping regression approach is used to decompose the relationships between the vegetation sensitivity (NDVI–PDSI regression slopes) and the principle climate factors (temperature and precipitation) associated with the drought. Significantly increased sensitivity to drought in warmer locations indicates the important role of temperature in exacerbating vulnerability; however, spatial precipitation variations do not demonstrate significant effects in modulating drought sensitivity. Based on annual NDVI response, chaparral is the most vulnerable community to warming, which will probably be severely affected by hotter droughts in the future. Drought sensitivity of coastal sage scrub (CSS) is also shown to be very responsive to warming in fall and winter. Grassland and developed land will likely be less affected by this warming. The sensitivity of the overall vegetation to temperature increases is particularly concerning, as it is the variable that has had the strongest secular trend in recent decades, which is expected to continue or strengthen in the future. Increased temperatures will probably alter vegetation distribution, as well as possibly increase annual grassland cover, and decrease the extent and ecological services provided by perennial woody Mediterranean climate ecosystems as well. divide this domain into the lower coastal plains ( < 400 m a.s.l.) in southwest and the relatively higher inland desert (600–1200 m a.s.l.) in northeast. The lowest elevations ( < 200 m a.s.l.) are dominated by developed or cultivated lands. The native Southern California forest, shrubland, and grassland are mainly distributed over the mountain ranges.


Introduction
Globally, CMIP5 climate models consistently project future warming and drying in most Mediterranean climate regions through the poleward-expanding subtropical subsidence [1][2][3]. Previous studies have suggested that the impacts of climate fluctuations on ecosystem processes become less predictable as climate extremes increase in frequency [4,5]. The resulting ecosystem changes can have many implications. As vegetation distributions change, so will habitat characteristics, potentially to the detriment of some species. Mediterranean climate regions are diversity hotspots, with large numbers of The land cover distribution reflects the rugged topography and the diversity in climate, with the more forested area at the tops of the Transverse Ranges; principally developed, cultivated lands or oak woodlands, shrublands and grasslands over the lower plains; and vast desert scrubs over the Mojave Desert ( Figure 2). The coastal mountains are dominated by CSS and chaparral communities, both of which are mainly composed of drought-tolerant species. However, CSS plants are generally shallow-rooted and drought-deciduous, compared to the mostly evergreen chaparral [39]. Due to the extremely dry summers, Southern California is a fire-prone region. Fire is a significant vegetation disturbance, which can degrade chaparral communities to grasslands abruptly. Shrub canopies recover in the following years; nevertheless, the post-fire regeneration largely depends on the fire intensity [57,58]. The land cover distribution reflects the rugged topography and the diversity in climate, with the more forested area at the tops of the Transverse Ranges; principally developed, cultivated lands or oak woodlands, shrublands and grasslands over the lower plains; and vast desert scrubs over the Mojave Desert ( Figure 2). The coastal mountains are dominated by CSS and chaparral communities, both of which are mainly composed of drought-tolerant species. However, CSS plants are generally shallow-rooted and drought-deciduous, compared to the mostly evergreen chaparral [39]. Due to the extremely dry summers, Southern California is a fire-prone region. Fire is a significant vegetation disturbance, which can degrade chaparral communities to grasslands abruptly. Shrub canopies recover in the following years; nevertheless, the post-fire regeneration largely depends on the fire intensity [57,58].

Vegetation Index Dataset
Satellite-derived Normalized Difference Vegetation Index (NDVI) data at a spatial resolution of 250 m and a temporal resolution of 16 days were acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, distributed by NASA's Earth Observing System Data and Information System (https://reverb.echo.nasa.gov/). The time span of the MOD13Q1 NDVI products (Collection 6, tile h08v05) cover from 18 February 2000 to 22 April 2017 [59]. The data after spring 2017 were excluded from this study, because of the extreme fire seasons in 2017 and 2018 (e.g., the Thomas Fire [60] and the Woolsey Fire [61]). These fires abruptly decreased the NDVI of large areas and hindered the analysis of climatic and drought effects.
We conducted some preprocessing for the NDVI data to ensure the reliability of this study. The pixels with negative NDVI values, which denote little or no vegetation (e.g., completely built-up urban areas), were removed before analysis. As snow tends to lower the NDVI [62], MODSCAG snow covered-area data [63], provided by NASA Jet Propulsion Laboratory, were applied to mask the NDVI pixels that had average snow fractions of more than 1% during a 16-day period. To match the temporal resolution of the drought index, the 16-day NDVI time-series were resampled to monthly time-series, by averaging the available tiles in each calendar month at each pixel.
Parameter-elevation Regression on Independent Slopes Model (PRISM, http://prism.oregonstate.edu; [68]) climate data serves as the foundational data set for computing the WWDT drought index. United States Historical Climatology Network (USHCN) data and Penn State STATSGO soil water holding capacity data also contributed to the climate data adjustment and the soil water balance calculations respectively [64]. The WWDT drought data has been intensively applied to climatological and hydrological studies in recent years (see [69][70][71][72][73][74][75]), and has been proved with good reliability. To match the spatial resolution of the MODIS NDVI maps, the PDSI dataset was interpolated to a 250 m × 250 m resolution using bilinear interpolation. The uncertainties and limitations of the data rescaling are discussed in Section 4.

Vegetation Index Dataset
Satellite-derived Normalized Difference Vegetation Index (NDVI) data at a spatial resolution of 250 m and a temporal resolution of 16 days were acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, distributed by NASA's Earth Observing System Data and Information System (https://reverb.echo.nasa.gov/). The time span of the MOD13Q1 NDVI products (Collection 6, tile h08v05) cover from 18 February 2000 to 22 April 2017 [59]. The data after spring 2017 were excluded from this study, because of the extreme fire seasons in 2017 and 2018 (e.g., the Thomas Fire [60] and the Woolsey Fire [61]). These fires abruptly decreased the NDVI of large areas and hindered the analysis of climatic and drought effects.
We conducted some preprocessing for the NDVI data to ensure the reliability of this study. The pixels with negative NDVI values, which denote little or no vegetation (e.g., completely built-up urban areas), were removed before analysis. As snow tends to lower the NDVI [62], MODSCAG snow covered-area data [63], provided by NASA Jet Propulsion Laboratory, were applied to mask the NDVI pixels that had average snow fractions of more than 1% during a 16-day period. To match the temporal resolution of the drought index, the 16-day NDVI time-series were resampled to monthly time-series, by averaging the available tiles in each calendar month at each pixel.

Drought Index and Climate Variables
The monthly gridded, self-calibrated Palmer Drought Severity Index (PDSI) data set with a spatial resolution of 4 km was downloaded from the WestWide Drought Tracker website ( [64]; WWDT, https://wrcc.dri.edu/wwdt/). PDSI is the most prominent index of meteorological drought in the United States [65], which incorporates antecedent precipitation, moisture supply, and moisture demand into a hydrologic accounting system [66]. Negative (positive) PDSI values indicate a drought (wet) period, and PDSI <−2, <−3, and <−4 refer to moderate, severe, and extreme drought, respectively [67]. Parameter-elevation Regression on Independent Slopes Model (PRISM, http://prism.oregonstate. edu; [68]) climate data serves as the foundational data set for computing the WWDT drought index. United States Historical Climatology Network (USHCN) data and Penn State STATSGO soil water holding capacity data also contributed to the climate data adjustment and the soil water balance calculations respectively [64]. The WWDT drought data has been intensively applied to climatological and hydrological studies in recent years (see [69][70][71][72][73][74][75]), and has been proved with good reliability. To match the spatial resolution of the MODIS NDVI maps, the PDSI dataset was interpolated to a 250 m × 250 m resolution using bilinear interpolation. The uncertainties and limitations of the data rescaling are discussed in Section 4.
The 30-year monthly PRISM Normals of temperature and precipitation were used in this study to demonstrate the normal climate conditions over the space of the study area [55]. This data set covers Remote Sens. 2019, 11, 2902 6 of 23 the period of 1981-2010 and has an initial spatial resolution of 800 m, which was rescaled on 250 m grids using bilinear interpolation. In addition, the monthly Daymet climate data set (version 3) on a 1 km grid was used to elucidate the climate change background for the period 1980-2017 in Southern California ( [76]; https://daymet.ornl.gov/). The long-term trends of monthly averages of maximum temperature, minimum temperature, as well as monthly total precipitation were analyzed.

Ancillary GIS Data
Fire perimeter data for the period January 2000-April 2017 were provided by the Department of Forestry and Fire Protection of California's Fire and Resource Assessment Program (FRAP, http: //frap.fire.ca.gov/). The FRAP compiled this statewide spatial database of fire perimeters from BLM (US Bureau of Land Management), NPS (National Park Service), and USFS (US Forest Service) fires. The FRAP fire database represents the best digital record of fire history in California, and its data quality is higher than that of MODIS fire data or other data sets [77,78]. The California State Wildlife Action Plan (SWAP; https://www.wildlife.ca.gov/SWAP) data were applied to identify the vegetation polygons of chaparral, CSS, grassland, forest, and developed land over Southern California [79]. The SWAP land cover data provides information on the distribution of vegetation types with high spatial resolution (30 m resolution) ( Figure 2). We follow [80] and the SWAP data were aggregated to a 250 m × 250 m resolution using majority statistics. This means that the vegetation type for each MODIS pixel was identified as the type that had the highest frequency at 30 m resolution in each 250 m × 250 m grid. This approach better reflects the regional vegetation distribution than the nearest-neighbor interpolation, and is less affected by data noise and misclassifications [80]. Additionally, a rigorous data quality control was implemented by the California Department of Fish and Wildlife in generating the SWAP database and, thus, it is the "best available" land cover data for California at present [50]. Vegetation cover typical of chaparral is comprised of sclerophyllous, evergreen-leaved shrubs with hard stems and deep roots with a variety of species. CSS is dominated by softer-leaved, woody shrubs and sub-shrubs. The rooting zone is shallower than chaparral and the dominant shrubs may be summer deciduous or possess smaller leaves during the dry season. Grasslands are typically dominated by a variety of non-native annual grasses, which die in the summer months and the cover regenerates through seed germination in the winter. Lower-elevation stands of woodland-forests are dominated by evergreen and deciduous oak, while the higher elevations (>1300 m a.s.l.) are mainly occupied by conifer stands (e.g., pines and firs).

Methodology
The monthly time-series of NDVI anomalies at the 250 m pixel level of resolution was calculated for the period from February 2000 to April 2017. Then, Pearson's correlations and linear regressions were performed between the NDVI anomalies and PDSI for each pixel. The NDVI-PDSI regression slopes represent the vegetation sensitivity to drought, and the NDVI-PDSI correlation coefficients indicate the strength of the relationships between the vegetation greenness and water availability. Correlation and regression analyses between NDVI and PDSI were also calculated for different seasons: winter, December-February; spring, March-May; summer, June-August; and fall, September-November. All the correlation and regression coefficients were tested at a 95% confidence level.
To detect long-term trends in the browning or greening of the vegetation over the past 17 years, a Theil-Sen estimator was utilized to measure the changing trends of the pixel-wise NDVI. The Theil-Sen algorithm is a robust trend detection approach, which gives the median changing slope of a time-series and, thus, is less affected by outliers than parametric approaches [81,82]. Though the NDVI time-series also contains some fluctuations induced by short-term disturbances (e.g., fire or diseases), the Theil-Sen method is appropriate for capturing the long-term greenness trends attributed to climate change. The Mann-Kendall test was applied to identify the significance of the NDVI trends. As a non-parametric test, the Mann-Kendall test is a reasonable method for detecting monotonic trends in series of environmental data, climate data, or hydrological data, regardless of the distribution of the data [83][84][85]. Given that autocorrelations commonly exist in environmental time-series, we used the pre-whitening algorithm proposed by [86] to remove the lag-1 serial correlations in the anomaly NDVI time-series before trend detection. Existence of a serial correlation can increase the probability that the Mann-Kendall test detects a significant trend, leading to a disproportionate rejection of the null hypothesis of no trend [87]. We also applied the same trend detection procedure to the monthly time-series of maximum temperature, minimum temperature, and precipitation, which indicated the recent climate change background in Southern California.
To compare the relative sensitivity of the five vegetation communities to drought, a dummy variable regression was performed, which uses the regression slopes between NDVI anomalies and PDSI as the dependent variable and vegetation as the unique predictor. Vegetation in developed land was selected as the reference category, as it showed the smallest seasonal variations. The dummy regression model gives the relative changes in NDVI-PDSI regression slopes for forest, chaparral, CSS, and grassland, as compared to developed land. By comparing the drought sensitivity of natural and developed-land vegetation, we may estimate how human activities mitigate the drought impact on urban vegetation and verify the assumption that natural vegetation is more vulnerable to drought than vegetation in developed lands.
Moreover, we employed multiple linear regression to further investigate the geographically different vegetation response to drought caused by climate, using the regression slopes between NDVI anomalies and PDSI as the dependent variable, while normal temperature and precipitation were the independent variables. As multiple linear regression assumes the independent variables are not highly correlated with each other, we used variance inflation factor (VIF) to detect multicollinearity [88]. A general rule is that the VIF should not exceed 10 to avoid severe multicollinearity [89]. Based on the two-fold regression, we can identify the impact of changing temperature and precipitation on the drought sensitivity of the vegetation types. To remove the influence of wildfire, we discarded all burned pixels over the study period, which were identified according to the FRAP fire perimeter data. As there are about 1 million MODIS pixels in this area, the large sample size may result in overestimation of minuscule effects in statistical analysis [90], represented by extremely small p-values. To overcome this problem, two solutions were utilized: first, a bootstrapping approach was applied to conduct the multiple linear regression, and within each of the total 1000 replicates, only 500 pixels were randomly selected to derive the regression coefficients and p-values; second, the standardized regression coefficient (slope) for each predictor in each bootstrapping replicate was also calculated as an effect size index [91]. The standardized slopes were calculated as where β j are the raw regression slopes and SD(X j ) and SD(Y) are the standard deviations of the independent variables (normal temperature and precipitation) and dependent variables (regression slopes between NDVI anomalies and PDSI), respectively. Standardized slopes in multiple regression analysis are scale-free estimates of the effect of a focal predictor on a target outcome [92]. These estimates are comparable in magnitude within different models and studies, and thus have advantages over partial correlation coefficients [93]. Therefore, the standardized slopes are more appropriate for demonstrating the impacts of the climate gradient on vegetation drought-tolerance than the raw regression slopes, and facilitate the comparison between temperature and precipitation contributions.

Vegetation Dynamics and Relationships with Drought
As indicated by the decreasing PDSI, the Californian climate has become drier over the past 30 years, where the relative drying on the south coast has been stronger than the rest of the state ( Figure 3). An extreme water deficit can exert negative impacts on ecosystem health, as the observed decline in vegetation greenness. Due to the persistent drying ( Figure 3 regions showed the most abrupt NDVI declines at the rate more than 0.009/yr ( Figure 4). However, browning was not constrained to burned sites. Besides the fire disturbances, the severe browning along the Californian south coast could mainly be attributed to climate disturbances. The spring (MAM) and summer (JJA) images show the greatest NDVI declines. The most severely browned regions include west Riverside County (more than 0.009/yr in spring), southwest Kern County (more than 0.007/yr in spring), south Los Angeles County (0.003-0.005/yr in spring and summer), and the western Transverse Ranges (more than 0.005/yr in all seasons).
Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 24 in vegetation greenness. Due to the persistent drying ( Figure 3) and frequent wildfires ( Figure 1) in Southern California, the natural vegetation might be expected to display general browning. Pixelwise Theil-Sen's slopes of NDVI represent the greening or browning of the vegetation ( Figure 4). The burned regions showed the most abrupt NDVI declines at the rate more than 0.009/yr ( Figure 4). However, browning was not constrained to burned sites. Besides the fire disturbances, the severe browning along the Californian south coast could mainly be attributed to climate disturbances. The spring (MAM) and summer (JJA) images show the greatest NDVI declines. The most severely browned regions include west Riverside County (more than 0.009/yr in spring), southwest Kern County (more than 0.007/yr in spring), south Los Angeles County (0.003-0.005/yr in spring and summer), and the western Transverse Ranges (more than 0.005/yr in all seasons). Though all the five major vegetation types showed stronger NDVI decline (Table 1), there were some seasonal differences between them. CSS (-0.0083/yr) and grassland (-0.0078/yr) showed greater NDVI declines in spring than the other seasons and vegetation communities. Chaparral showed strong browning across all four seasons. The fact that chaparral and forest did not display apparent seasonal variations in the NDVI trends might be related to their deep-rooted features. CSS and grassland communities displayed less browning than chaparral and forest in summer and fall. The chaparral communities (-0.0078--0.0082/yr) and the developed-land vegetation (-0.0033--0.0027/yr) showed the largest and smallest NDVI declines, respectively, throughout the year (see Table 1).   The significant correlations and regressions between NDVI anomalies and PDSI for all seasons in nearly all pixels suggest that the vegetation was very sensitive to drought in this region ( Figures 5  and 6). The strongest correlations were found in the summer months (June-August), while the highest regression slopes were shown in winter (December-February). The Mediterranean climate Though all the five major vegetation types showed stronger NDVI decline (Table 1), there were some seasonal differences between them. CSS (−0.0083/yr) and grassland (−0.0078/yr) showed greater NDVI declines in spring than the other seasons and vegetation communities. Chaparral showed strong browning across all four seasons. The fact that chaparral and forest did not display apparent seasonal variations in the NDVI trends might be related to their deep-rooted features. CSS and grassland communities displayed less browning than chaparral and forest in summer and fall. The chaparral communities (−0.0078-−0.0082/yr) and the developed-land vegetation (−0.0033-−0.0027/yr) showed the largest and smallest NDVI declines, respectively, throughout the year (see Table 1). The significant correlations and regressions between NDVI anomalies and PDSI for all seasons in nearly all pixels suggest that the vegetation was very sensitive to drought in this region ( Figures 5  and 6). The strongest correlations were found in the summer months (June-August), while the highest regression slopes were shown in winter (December-February). The Mediterranean climate can account for this seasonal pattern (i.e., most of the precipitation and vegetation growth occurs in winter, and all the natural plants must endure water scarcity in summer). In consequence, moderate magnitudes of the correlation (0.47-0.62) and regression coefficients (0.0109-0.0204) between NDVI and PDSI for the natural vegetation were observed for spring (March-May) and fall (September-November) as transitional periods (Tables 2 and 3).   Dark gray lines refer to burned areas. Insignificant pixels have been excluded. White areas indicate data gaps due to snow and cloud. Correlation coefficients are the largest in summer and smallest in winter.   Although desert vegetation is adapted to low moisture availability, the vegetation greenness in the Mojave Desert (northeast study area) showed surprisingly strong positive correlations with PDSI, indicating vegetation sensitivity to moisture stress there ( Figure 5). However, the regression slopes were smaller compared to the coast, perhaps due to lower relative sensitivity and the already scarce vegetation and lower NDVI value in the desert ( Figure 6). Weak correlations and high regression coefficients were found in the burned regions, which can be attributed to abrupt fire disturbances (Figures 5 and 6). For example, wildfires in drought years may largely reduce the NDVI of the burned areas, and when precipitation rebounds after drought, post-fire regeneration can bring a rapid increase in NDVI. Besides the desert and burned areas, the vegetation in the coastal mountain ranges demonstrated striking relations with the PDSI (e.g., in the San Rafael Mountains, Santa Ynez Mountains, Santa Monica Mountains, Santa Ana Mountains, and Riverside-Perris Plain).
We found that most of the browning regions showed close relationships between NDVI anomalies and PDSI (Figures 4 and 5). This is consistent with the general browning in Southern California due to decreasing PDSI during recent years ( Figure 3). As a measurement of dryness based on precipitation and temperature, reduction of precipitation and elevated temperature can both decrease PDSI. A large portion of Southern California has shown an appreciable warming trend in all seasons over the past 37 years (Figures 7 and 8). A few regions near the San Emigdio Mountains displayed cooling in the night (decreased minimum temperature), which was offset by stronger daytime warming. At the same time, precipitation did not generally demonstrate significant changes, with only the southwest Mojave Desert showing decreased precipitation in spring (Figure 9). We found that most of the browning regions showed close relationships between NDVI anomalies and PDSI (Figures 4 and 5). This is consistent with the general browning in Southern California due to decreasing PDSI during recent years ( Figure 3). As a measurement of dryness based on precipitation and temperature, reduction of precipitation and elevated temperature can both decrease PDSI. A large portion of Southern California has shown an appreciable warming trend in all seasons over the past 37 years (Figures 7 and 8). A few regions near the San Emigdio Mountains displayed cooling in the night (decreased minimum temperature), which was offset by stronger daytime warming. At the same time, precipitation did not generally demonstrate significant changes, with only the southwest Mojave Desert showing decreased precipitation in spring (Figure 9).

Relative Sensitivity of the Major Vegetation Types
Due to varied responses of the five vegetation communities to drought and increasing temperatures/aridity, the regression slopes between NDVI anomalies and PDSI varied among vegetation types. The relative sensitivity of the vegetation also displayed seasonal differences (Figure

Relative Sensitivity of the Major Vegetation Types
Due to varied responses of the five vegetation communities to drought and increasing temperatures/aridity, the regression slopes between NDVI anomalies and PDSI varied among vegetation types. The relative sensitivity of the vegetation also displayed seasonal differences ( Figure  10). Compared to developed land, grassland and CSS showed significantly increased drought

Relative Sensitivity of the Major Vegetation Types
Due to varied responses of the five vegetation communities to drought and increasing temperatures/aridity, the regression slopes between NDVI anomalies and PDSI varied among vegetation types. The relative sensitivity of the vegetation also displayed seasonal differences ( Figure 10). Compared to developed land, grassland and CSS showed significantly increased drought sensitivity in winter and spring, and a mild increase of drought sensitivity in summer and fall ( Figure 10 and Table 3). Grassland had the largest seasonal fluctuations in response to drought stress (Table 3). Chaparral showed the highest drought sensitivity in the warm seasons (median p < 0.05, Figure 10 and Table 3). Forest displayed an insignificant increase in drought sensitivity, compared to developed land, in the cold seasons, but a significantly increased sensitivity from summer to fall ( Figure 10).

Climatic Effects on Shifting the Drought Impacts
The spatial variations of the NDVI-PDSI regression slopes of the four natural vegetation types against the normal temperature and precipitation gradients suggest that all the regression slopes increased in geographic regions with higher temperatures but remained stable across geographic precipitation gradients (Figure 11). In other words, natural vegetation in the warmer regions of the study area were more sensitive to drought, as indicated by their larger NDVI-PDSI regression slopes. However, the spatial precipitation patterns did not demonstrate a similar geographic influence. Additionally, the vegetation in developed lands seems to be less affected by temperature or precipitation changes ( Figure 11). The seasonal correlation coefficients between NDVI anomalies and PDSI reflected the general water stress variations of the five vegetation types throughout the year (Table 2). Winter (summer) showed the lowest (highest) median correlation coefficients, indicating less severe (extreme) drought stress due to favorable (limited) water conditions. Spring and fall displayed intermediate correlation coefficients, suggesting moderate drought stress. Moreover, the NDVI-PDSI correlations also displayed slight differences among the five ecosystems (Table 2). Chaparral and CSS showed relatively higher correlation coefficients in all seasons, indicating that the two communities had a more consistent response to drought in time, as compared to forest, grass, and developed land.

Climatic Effects on Shifting the Drought Impacts
The spatial variations of the NDVI-PDSI regression slopes of the four natural vegetation types against the normal temperature and precipitation gradients suggest that all the regression slopes increased in geographic regions with higher temperatures but remained stable across geographic precipitation gradients (Figure 11). In other words, natural vegetation in the warmer regions of the study area were more sensitive to drought, as indicated by their larger NDVI-PDSI regression slopes. However, the spatial precipitation patterns did not demonstrate a similar geographic influence. Additionally, the vegetation in developed lands seems to be less affected by temperature or precipitation changes ( Figure 11). Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 24 Figure 11. Quartile statistics of the regression slopes between NDVI anomalies and PDSI along the gradients of normal temperature (left) and precipitation (right) for chaparral (1st panel), coastal sage scrub (2nd panel), grassland (3rd panel), forest (4th panel), and developed land (5th panel). NDVI-PDSI regression slopes increase with temperature gradients, but are stable as precipitation changes.
The bootstrapped regressions indicate that the standardized slopes (beta) against precipitation were smaller than those of temperature for all the major vegetation types (Figure 12). We calculated the variance inflation factor (VIF) for all the multiple linear regressions ( Figure S1). All the VIFs were <1.5, and thus these regressions did not violate the assumption of no multicollinearity [89]. The beta values of precipitation mainly ranged between -0.3 and 0.3, which denotes that one SD (standard deviation) precipitation change would only induce <0.3 SD change of the regression slopes between NDVI anomalies and PDSI. In contrast, a 1 SD temperature increase would raise the NDVI-PDSI regression slopes by up to 0.5 SD. Additionally, precipitation showed opposite effects in changing the vegetation sensitivity between the wet and dry seasons. In other words, more precipitation could Figure 11. Quartile statistics of the regression slopes between NDVI anomalies and PDSI along the gradients of normal temperature (left) and precipitation (right) for chaparral (1st panel), coastal sage scrub (2nd panel), grassland (3rd panel), forest (4th panel), and developed land (5th panel). NDVI-PDSI regression slopes increase with temperature gradients, but are stable as precipitation changes.
The bootstrapped regressions indicate that the standardized slopes (beta) against precipitation were smaller than those of temperature for all the major vegetation types (Figure 12). We calculated the variance inflation factor (VIF) for all the multiple linear regressions ( Figure S1). All the VIFs were <1.5, and thus these regressions did not violate the assumption of no multicollinearity [89]. The beta values of precipitation mainly ranged between −0.3 and 0.3, which denotes that one SD (standard deviation) precipitation change would only induce <0.3 SD change of the regression slopes between NDVI anomalies and PDSI. In contrast, a 1 SD temperature increase would raise the NDVI-PDSI regression slopes by up to 0.5 SD. Additionally, precipitation showed opposite effects in changing the vegetation sensitivity between the wet and dry seasons. In other words, more precipitation could decrease the NDVI-PDSI slopes in winter and spring, and increase them in summer and fall, resulting in insignificant annual effects (p > 0.05) for most vegetation communities, except for grassland ( Figure 12). Temperature demonstrated significant (p < 0.05) and positive beta values in nearly all bootstrapped regressions (seasonal and annual), except for grassland in spring and summer, as well as developed land in all seasons. Thus, warmer temperatures consistently increased the natural vegetation sensitivity to drought, and the effects were prone to be greater in the dry season. Strikingly, the chaparral communities displayed very high median standardized slopes of temperature (>0.3 SD, p < 0.05) in all seasons (Figure 12e). decrease the NDVI-PDSI slopes in winter and spring, and increase them in summer and fall, resulting in insignificant annual effects (p > 0.05) for most vegetation communities, except for grassland ( Figure  12). Temperature demonstrated significant (p < 0.05) and positive beta values in nearly all bootstrapped regressions (seasonal and annual), except for grassland in spring and summer, as well as developed land in all seasons. Thus, warmer temperatures consistently increased the natural vegetation sensitivity to drought, and the effects were prone to be greater in the dry season. Strikingly, the chaparral communities displayed very high median standardized slopes of temperature (>0.3 SD, p < 0.05) in all seasons (Figure 12e).

Discussion
Our analysis reveals the relative sensitivity of five major wildland vegetation types and developed areas to increased warming and aridity in Southern California, a Mediterranean climate region. Climatic water deficit (CWD), which combines the influences of temperature and precipitation, has been shown to be an important driver of woody plant decline and mortality in drought [94], which is consistent with the results presented here. Quetin and Swann [95] recently investigated the impacts of changing temperature and precipitation on global NDVI. Here, we further focus on how a changing climate will affect vegetation drought sensitivity. The bootstrapped multiple regression model utilized in our study dissociated the effects of temperature and precipitation on drought stress, and allowed for examination of the seasonal variations of these effects. The bootstrapping approach successfully reduced the overestimation of minuscule effects, which generally occurs in big data (e.g., remotely sensed data) analytics [90]. Precipitation consistently showed more insignificant results in predicting vegetation sensitivity to drought, while temperature displayed significant effects in nearly all bootstrapped regressions.
Previously, [23] has found that the forests in southwest U.S. are particularly sensitive to the increasing aridity and warmth. We further compared the differential responses of drought sensitivity (NDVI-PDSI regression slopes) to temperature and precipitation changes in forest, shrubs, grass, and developed land. All the natural vegetation types showed greater sensitivity in the warmer areas of the region, implying that temperature increases alone will be an important determinant of vegetation vulnerability to global warming in the 21st century. Vegetation in developed lands were consistently less responsive to warming, indicating the effects of irrigation in droughts, although winter precipitation increases seem to be able to reduce drought sensitivity in urban vegetation (Figure 12f). Grasslands showed less sensitivity in response to the combined effects of high temperatures and drought stress; while chaparral displayed very high sensitivity in all seasons. Thus, grassland and chaparral seem to be the most resistant and vulnerable communities, respectively, to drought under the effects of a warmer climate (Figure 12e). This is consistent with a recent study by [96], who concluded that forest and grass displayed longer and shorter drought legacy effects, respectively, but both showed greater drought resilience than shrubs. These differences arise from different eco-hydrological and physiological responses of these plant's functional groups (e.g., rooting systems, water-use patterns, and differential hydraulic responses).
It should be mentioned that uncertainties exist in the relationship between extreme events in climate and vegetation (e.g., non-linearities and time lags), but pixels having more drought stress should have experienced more vegetation extremes and, therefore, the slopes of the NDVI-PDSI linear regressions could be used as an indicator of apparent sensitivity [51]. Another limitation of this study is related to the inconsistency of the databases used, in terms of spatial resolution. Bilinear interpolation was used in rescaling the PDSI and PRISM climate normals to MODIS resolution. For these continuous variables, bilinear interpolation can provide more accurate predictions than the other simple resampling approaches, such as nearest neighbor interpolation, and it has been widely used to rescale climate variables in similar studies, e.g., [97][98][99][100]. Some researchers have developed statistical or dynamical downscaling techniques for climate variables, e.g., [101,102]. However, as the space-time variability and patterns of a climate variable are highly dependent on the location of this variable, it is not practical to directly apply a downscaling approach developed for another region to our study region, or to use a downscaling approach for all climate variables. In addition, as PDSI reflects the variations of many climate variables (e.g., temperature, precipitation, evaporation, soil moisture, and wind speed), there is still no ready-to-use downscaling method for PDSI. Moreover, [17] analyzed how spatial resolution changes affect the relationship between PDSI and NDVI, and they found that the NDVI-PDSI correlations and regression coefficients did not change much when the data resolution was adjusted from 4 km to 500 m. Besides the spatial resolution issue, uncertainties may arise from the differential temporal resolution of PDSI and NDVI data. To match the PDSI time step, we adjusted the MODIS NDVI data set from 16-day to monthly. The higher frequency of the NDVI observations can capture more vegetation dynamics. However, as this study was focused on climatological analysis, and as the three droughts in Southern California during 2000-2017 were all prolonged droughts (2-4 years), the degraded temporal resolution of MODIS NDVI should only have little influence on the results included in this paper. Thus, we assume the uncertainties and biases induced by the different data resolutions cannot largely change the conclusions presented here. In addition, the lagged relationship between hydrological variables (precipitation, soil moisture, and drought) and vegetation growth also brings uncertainties to this study. Compared to the newly developed drought indices, such as the standardized precipitation index (SPI) and the standardized precipitation evapotranspiration index (SPEI), PDSI is limited by its fixed time scale of one month. However, the PDSI of a given month can represent an accumulation of drought over the previous few months and is highly correlated with the top 1 m soil moisture [103,104]. SPI is calculated solely based on precipitation and does not take temperature influence into account. SPEI is better than PDSI in reflecting the variations of soil moisture and evapotranspiration in drought and is superior for its varying time scales [103]. As a traditional drought index, PDSI has long been used for drought measurement and has advantages in providing comparable results with other related studies, e.g., [15,16,105,106]. As [107] suggested, no single drought index can capture all aspects of soil moisture drought. An interesting work in the future would be investigating the vegetation responses to drought and climate change using other drought indices and varying time scales.
Previous studies (e.g., [108][109][110][111]) have suggested that the sprawling development pattern in Southern California has been an important driver of contemporary chaparral conversion to herbaceous cover, both through the direct removal and fragmentation of habitat, but also through its indirect role in driving annual grass expansion associated with shortened fire intervals. As we have excluded the burned areas in this study, the NDVI change analyses should not be affected by the wildfires. As the SWAP land cover map is a recently updated (in 2015) database, the areas newly occupied by residents were identified as developed lands, instead of as natural vegetation. Thus, the NDVI trends and vegetation sensitivity for the developed lands in this study may include some signals of type conversion from natural vegetation communities to urban vegetation. However, as our study is limited to the MODIS era (after 2000), the uncertainties due to vegetation type conversions are restrained and should not alter the main conclusions.
Relative differences in the drought and temperature sensitivities of the five wildland vegetation types and their related seasonal changes in NDVI also reflect their different phenological characteristics. For example, the annual Californian grassland and summer deciduous or dimorphic CSS are well-adapted to the summer drought of the Mediterranean climate; as they die back, they become dormant or reduce photosynthesis when the soil moisture significantly decreases in summer [40]. In the cool and humid winter, the two communities reach maximum biomass quickly. By contrast, chaparral and forest trees consist of many evergreen species, which carry on photosynthesizing over the summer months and show less seasonality [40]. As a result, the CSS and grassland flora displayed higher NDVI-PDSI regression slopes than the chaparral and forest communities in winter and spring and lower slopes from summer through fall, partly because of the seasonally low NDVI of the two former communities in the warm season ( Figure 10).
Climate trend detection suggests an intense warming in the region since 1980s, especially in summer and fall, when the most extreme seasonal water shortage occurs (Figures 7 and 8). At the same time, precipitation did not display any meaningful trends, except for a slight decrease over the Southern Mojave Desert in spring. Previous studies suggest that the differing drought intensities and inter-annual variations of precipitation in Northern and Southern California are related to a hydroclimatic dipole [17,112]. Both the historical observations and the future projections support a more severe drying in Southern California than the north [17,113]. Within Southern California, it is still not clear what factors are controlling the spatial patterns of the drying and warming trends. Both the large-scale atmospheric circulation and the regional-scale land surface conditions (e.g., topography, albedo, soil moisture, and land covers) are potential influences.
The recent warming seems to be an important component in any long-term trends in browning over Southern California (Figure 4). There were also a few regions at higher elevations (>1500 m a.s.l.) that demonstrated greening trends over the past 17 years (Figure 4). The greening of the cold-limited montane vegetation (mainly in winter) might be caused by climate warming (i.e., due to the favorable temperature and earlier snowmelt) [114]. In addition, it is still not clear why some highly elevated areas over the San Emigdio Mountains showed significant night cooling and daytime warming (Figure 8). As a fire-prone region (Figure 1), it should be mentioned that an abrupt vegetation removal by wildfires may create a more arid environment with enhanced diurnal temperature fluctuations, i.e., increasing the temperature in daytime but decreasing it at night [54]. Moreover, since chaparral shrubs are susceptible to freezing injury [115], persistent night cooling of the San Emigdio Mountains may be as harmful as high temperatures and drought to the wildely distributed chaparral plants in this region.
Jacobsen and Pratt [116] found that the Southern California shrublands are relatively resilient to short-term high-intensity droughts, but displayed significant dieback and mortality during the recent prolonged drought (2012-2016); they concluded that the highly cavitation-resistant shrubs are the ones most affected by long-term high-intensity drought. This is supported by [50], who argued that the deep-rooting system of chaparrals buffers them from summer drought but makes them particularly vulnerable to multi-year drought, where the soil water stored in the deep layers is reduced by persistent evaporation. Our results further reveal that the chaparral communities seem to be less resilient to climate disturbance than the other vegetation types. They will likely be impacted more severely by future climatic warming and drying, due to their highest overall sensitivity to drought stress exacerbated by annual temperature increases (Figure 12e). Previously, [50] found that the cooler and moister climate at high elevations of the Californian south coast could produce decreased sensitivity of chaparral to drought but did not change the sensitivity of CSS. However, coastal fog might weaken the temperature effects on drought stress in CSS [50]. In this paper, we extended the study area to the inland Southern California mountains, and high temperatures, indeed, displayed a significant effect of increasing drought impacts on CSS in all seasons ( Figure 12). Thus, the drought stress of the CSS plants, particularly where fog is limited, will probably also increase under a warmer future climate. In the following work, we are going to apply this analysis to other Mediterranean-climate regions exposed to drought, and disentangle the effects of snow, fire, groundwater level changes, and low temperature extremes for the mountainous areas. Moreover, it is important to further reveal the role of root depth in modulating the drought sensitivity of the Mediterranean-climate ecosystems.
In response to the increasing aridity and frequent fires, the Southern Californian ecosystems may shift over time, perhaps favoring grasslands and generally lower biomass covers [14,117]. As previously noted, any replacement of woody vegetation by grasslands or developed lands will alter ecosystem services, through loss of certain habitat types in affected locales and increased vulnerability of slopes to erosions and landslides. Additionally, since chaparral shrubs may have an important role in holding soil water and increasing groundwater recharge (provisioning) through channel infiltration, a degradation of the chaparral community may affect the local groundwater level [118,119]. In addition to these regional impacts, vegetation transitions could impact the carbon balance. Mediterranean climate vegetation communities have been shown to be among those displaying the highest net primary productivities (NPP) amongst the world biomes [120,121], but the net ecosystem carbon exchange of Mediterranean climate vegetation is highly sensitive to drought [122,123]. Consequently, a warmer and drier climate could significantly deteriorate the carbon budget over Mediterranean climate regions, such as Southern California. Decreased vegetation productivity coupled with the loss of woody vegetation would favor the conversion to grasslands and developed lands.

Conclusions
In this study, we utilized a bootstrapping regression to investigate the effects of changing temperature and precipitation in modulating the drought sensitivity of Mediterranean ecosystems in Southern California. The main conclusions are summarized as follows: (1) Temperature increase can significantly aggravate the susceptibility of the major Mediterraneanclimate vegetation types to drought. The recent greenness declines may be related to the synchronous compound influence of high temperatures and drought. (2) Precipitation changes display insignificant or less effects in adjusting the vegetation sensitivity to drought.
Chaparral seems to be the most vulnerable community to the future hot droughts in all seasons, and CSS will likely be more sensitive to drought from fall to winter under a warmer climate. Climate warming may only exert small effects in grassland and developed-land vegetation.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/24/2902/s1, Figure S1: Variance inflation factors (VIFs) of the multiple linear regressions shown in Figure 12. VIF is a widely used indicator for measuring the multicollinearity in multiple linear regressions. Normally, VIF>10 indicates severe multicollinearity, 5<VIF<10 indicates weak multicollinearity, and VIF<5 indicates no multicollinearity. Nearly all the VIFs in this study <1.5, indicating these multiple linear regressions are not affected by multicollinearity.