Analysis of the Ecological Effects of Decadal Large Scale Intermittent Annual Water Allocation using Satellite Data in Baiyangdian Wetland, Northern China

In this study, the ecological effects of intermittent water allocation with emphasis on spatiotemporal responses of the corresponding vegetation were analyzed using remote sensing data and GIS-based buffer technology considering the period from 1st July 2000 to 31st December 2009. Three sampling sites (Angzh, Wangk, and Xidayang) with different water flow paths and three buffer distances were distinguished in the research. The Seasonal-Trend decomposition procedure based on Regression (STR) trend extraction and its corresponding linear regression and anomaly detection were executed to determine temporal variations of vegetation under the effects of water allocation. ANOVA and PCA methods were employed to identify the spatial responses of vegetation to different water flow paths and buffer distances. The results were as follows: (1) NDVI except NDVImin displayed higher values during the period without water allocation; (2) extremely significant decline trends (p<0.001) of all NDVI categories were observed in all sites at all buffer distance levels, except for NDVImin at buffer distances of 2 km and 4 km in Angzh, showing stronger fluctuations of frequency after 2008 as well as the decline gradient with the extent of buffer distance to river. The anomaly detection results provided similar evidence of stronger NDVI fluctuations after 2008; (3) water allocation had extremely significant effects on regional vegetation coverage (p<0.01) with a decline gradient of statistical p values along enlarged buffer distances. Our results provide evidence of spatial and temporal differences in vegetation response to water availability due to the intermittent frequency water allocation implemented via different river channels. The findings of this study will deepen our understanding of the effects of water division on regional vegetation restoration and can be used to develop a practical strategy for effective implementation of water allocation.


Introduction
Vegetation cover is a critical component of wetland ecosystems in arid and semi-arid regions; however, water scarcity has remained a serious issue for a long time [1]. Vegetation dynamics are considered the primary ecological indicator of ecosystem response to water availability in wetland ecosystems. The dynamics of vegetation cover are dominated by the availability of water to a large extent in its temporal and spatial distribution [2]. Therefore, the relationship between vegetation and water availability is becoming increasingly important in the fields of global environmental change and ecological restoration.
Managing water for the sustainable development of ecosystems is both a technical and a governance challenge in which balancing water loss and vegetation restoration plays a central role [3][4][5]. Implementation of water allocation aimed at maintaining ecosystem sustainability is a mechanism that may be used in basin-scale water management [3]. In arid and semiarid watersheds, the contrast between high water usage in midstream areas (i.e. irrigation for agriculture) and water crisis in natural ecosystems along downstream areas has become increasingly prominent [6,7]. The higher consumption of water in the midstream area implies lower water availability for the downstream area, especially in areas with originally scarce water resources. The responses of wetland ecosystems to river flow regulation associated with water diversion involve many complex ecohydrological processes [8,9]. Therefore, it is crucial to determine an appropriate level of water diversion considering both human needs and ecological restoration. Understanding changes in vegetation before and after water diversion and adjusting the eco-hydrological linkages between midstream and downstream areas are essential for water resource planning in water scarce areas.
Water allocation has become a necessary measure to improve the reallocation of water resources for wetland ecosystem restoration in the Central Asia [10,11] and Northwestern China [12,13]. Remotely sensed data, such as normalized difference vegetation index (NDVI), serve as key indicators of vegetation status and are useful parameters in studies of terrestrial vegetation cover. The time series of NDVI data provide important information for evaluating vegetation change before and after water allocation to improve strategies of large-scale water division. The regional ecosystem of arid and semi-arid areas, including vegetation cover, is sensitive to fluctuations of local water resource availability [14][15][16]. The vegetation cover of wetland watersheds coupled with its eco-environmental quality and stability are crucial for maintaining a healthy ecosystem [17,18]. Understanding vegetation responses to combined water resource availability and the corresponding management planning would facilitate regional ecosystem management. Generally, ecosystem management strategies have been shown to have greater influence on vegetation distribution [19,20]. In the restoration of large-scale vegetation distributions, the effects of regional sustainable ecosystems could be reflected by NDVI [21,22]. Therefore, understanding gradients of NDVI variation and factors influencing NDVI response to water availability is important for sustainable ecosystem management. Omute et al. [23]revealed the strong interactions between NDVI and water levels especially in low precipitation years in Lake Victoria. Cao et al. [24] compared the water balance in potential natural vegetation restoration to guide future ecological restoration planning. Zhang et al. (2016a) emphasized on the importance of water allocation to vegetation restoration rather than conservation of natural vegetation in China. Previous research focused more on water allocation models considering the evaluation of ecosystem services [25], gaining public acceptance for larger water projects and policy changes [26], affecting the responses of vegetation dynamics to hydro-climatic factors [11], fulfilling water reclamation management [27], changing operation rules with hydrologic statedependent multi-reservoirs [3], as well as the relationships with soil, hydrology, vegetation and climate change [5,28,29]. These studies are mostly concentrated on the overall effects for a time period of water division or a specific effect during one water division, and on the final effects on the ecology and human society. Studies on the ecological effects of processes of intermittent water allocation, as well as the corresponding the strategy and spatial responses are lacking.
This study selected three sampling sites along water flow paths in the Baiyangdian watershed as examples of areas where constant water replenishment will be implemented in the future. These areas are the top priority for ecological conservation and restoration because the Xiongan New Area has been constructed as China's third state-level new area in April 2017. The sampling sites correspond to three water allocation paths with three upstream reservoirs; and they were selected considering less interference of river branches and ample distance from downstream lakes. With the anticipated population boom and the rapid economic development in the processes of regional construction, the consumption of water will increase dramatically and water availability for ecological processes will diminish. Determining the relative importance of NDVI categories at each site using PCA could enhance our understanding of water allocation strategy and the response of ecological processes to intermittent water division as well as spatial distance. This study used 10-year NDVI time series with 10-day syntheses datasets (1 st July 2000 to 31 st December 2009) of Baiyangdian wetland to examine the temporal and spatial response of NDVI to water resources. GIS-based three buffer distances (2 km, 4 km and 6 km distance to river channel center, respectively) were employed to assess the spatial response of NDVI to water resources. The NDVI series data were extracted to three distance gradient categories based on buffer levels at each site. The categories of each buffer distance level and each site included: the minimum value of NDVI of the all grid cells at the respective buffer distance level (NDVI min ); the mean value of NDVI (NDVI mean ); the maximum pixel value of NDVI, meaning the value of the largest number of all grid cells at the respective buffer distance level (NDVI most ); and the maximum value of NDVI (NDVI max ). ANOVA analysis was performed using the data for trends and anomaly changes in each NDVI category and then temporal and spatial vegetation responses to water allocation for periods before water allocation, during water allocation and after water allocation. Our main objectives were: (1) to elaborate decadal trends and abrupt patterns caused by intermittent water division in each NDVI category; (2) to identify temporal and spatial responses of vegetation to intermittent water allocation as well as the relative importance of NDVI categories; and (3) to propose a new perspective for implementation of water allocation at different time scales via different division paths.

Study area
The Baiyangdian watershed (38°10'-40°00' N and 113°40'-116°20' E), located 130 km south of Beijing (48°43'-39°02' N and 115°38'-116°07' E), is the largest natural wetland in the North China Plain (Fig. 1). In April 2017, Xiongan New Area as China's third state-level new area after Shenzhen Special Economic Zone and Shanghai Pudong New Area, achieved profound significance in economic and social development. Accordingly, the restoration of Baiyangdian wetland has become the top priority under the New District Construction. The wetland belongs to parts of the Haihe River basin, which is one of ten large basins in China. The total area of the Baiyangdian watershed is 31,200 km 2 and includes eight rivers: the Caohe, Fuhe, Juma, Pinghe, Puhe, Tanghe, Xiaoyi and Zhulong rivers. There are three reservoirs upstream of the Lake: Angezhuan reservoir (Angzh), Wangkuai reservoir (Wangk), and Xidayang reservoir (Xidayang). The topography does not vary greatly in such a flood plain. The regional climate is dominated by continental monsoons of humid or sub-humid region with an average rainfall of 530.4 mm, and precipitation mostly occurs in July and August. There is a distinct seasonality in the annual patterns of precipitation, evaporation, and air temperature. In the recent decades, the natural hydrologic regime of the upstream of the watershed reservoir, which is aimed at benefiting economic and social water users, has been profoundly changed. The annual natural outflow of the upstream reservoir is now inadequate to meet the environmental flow requirements of the downstream wetlands and their ecosystem health restoration. Water scarcity has become a crucial problem for the wetland restoration. The corresponding vegetation and agriculture have been greatly influenced. To alleviate water resource scarcity, water allocations from upstream reservoirs have been implemented since the 1990s (Table 1), resulting in significant impacts on regional hydro-ecological variations.

Data source
Series records of water allocation data for Baiyangdian watershed were obtained from the Anxin Environmental Bureau. The sources for the corresponding 10-year NDVI time series datasets (1 st July 2000 to 31 st December 2009) with a spatial resolution of 1 × 1 km were the PROBA-VGT S10 TOC NDVI products obtained from European Space Agency (http://www.vito-eodata.be). The VGT S10 products (10-day syntheses) were compiled by merging segments (data strips) acquired in a 10-day period. All segments of this period were compared pixel-by-pixel to select the 'best' surface reflectance values. These products provided data from all spectral bands, the NDVI, and auxiliary data on image acquisition parameters. The NDVI product ensured data quality with a minimum effect of cloud covers, and the pixel brightness count was the ground area's reflectance (corrected for atmospheric effects). All NDVI images were reprojected in UTM coordinates as metadata for the research.

Spatial buffer distance levels and vegetation ecological effects
The source of the corresponding water allocation records for decadal NDVI time series dataset (10-day syntheses with a spatial resolution of 1 km 2 ) were employed to determine the ecological effects of water allocation on vegetation. The corresponding flow paths of water divisions from upstream reservoirs was distinguished, as shown in Figure 1. Detailed information of water transformation from three upstream reservoirs is summarized in Table 1. To evaluate the magnitudes of ecological improvement attributable to the frequency of water allocation, three buffer distance levels (distances of 2 km, 4 km and 6 km from the river center, abbreviated as bf2, bf4, and bf6 respectively) were considered in the research. Areas in the basin with relatively few river branches and weak anthropogenic activities, as well as relatively gentle slope were considered during the selection of sampling sites (Angzh, Wangk, and Xidayang). Grid numbers with 1 km 2 spatial resolution of each grid at different buffer distances in the three sampling sites were summarized. The spatial buffer distance of each site and

STR and Trend component extraction
The seasonal variation in time series data is the most important characteristic and the basis of all seasonal adjustment procedures when building models to fit observation data. Various methods are available for handling seasonal components, such as Seasonal ARIMA (SARIMA), X-11-ARIMA [30] and X-12-ARIMA [31] and X-13ARIMA-SEATS [32], as well as STL (Seasonal-Trend decomposition using Loess), which has been widely used because of its availability in R [33]. However, only few of them afford the clarity, simplicity, and generality required for handling several problems that arise with seasonal data decomposition [34]. The STR (a Seasonal-Trend decomposition procedure based on Regression) could solves such problems by performing re-cast in the framework of ordinary least squares or quantile regression. In particular, the robust version of STR assumes a double exponential distribution for the residuals, and trend, seasonal and predictor coefficient changes compared to the normal distribution assumed in the STR model, thus providing detailed changes of the trend component.
In the simplest STR model, time series data may be considered the sum of three components: one highfrequency seasonal component, one low-frequency longterm component (or trend), and a residual component (variation not explained by time), which are expressed as: where T t is the trend, S sn(t), t is the additive seasonal component, which is a k×n matrix (k is number of seasons and n is length of the time series), and R t is the "remainder" component. The STR model uses the maximum likelihood to estimate cross-validated residuals of the seasonal components, and cross validation to estimate the residuals of minimizing SSE (the sum of squares due to error) by finding optimal parameters using R core "stats" package. Considering that the Robust STR can tolerate outliers well using the quantile regression approach (only 0.5 quantile is used), which could provide details of the trend component, we selected Robust STR to extract the trend component in the research. All STR analyses were performed using the "stR" package in R language (http://cran.R-project.org/). Because the data have 10-day granularity, the yearly seasonality included a period of 36 observations. Thus, gaps of 36 observations were used for cross validation, and 95% confidence intervals were calculated in the STR analysis. A more complete description of this methodology can be found in Dokumentov and Hyndman [34].

Linear regression and anomaly detection
The extracted trend components were measured as continuous time series variables, and the linear regression analysis was performed to explore possible statistical trends of the magnitude of NDVI changes over time. Specifically, the simple linear regression analysis was employed to explain changes in the different NDVI categories over time at each buffer distance level at the three sites. The partial least squares (PLS) regression was applied to capture the overall sensitivity of temporal variables in different years. The overall goodness-of-fit of the regression model was evaluated by the regression coefficients and the p-value, and the proportion of the observed variations in the response variables could be explained by the PLS model. The linear regression analyses were performed in R (http://cran.R-project.org/).
Anomaly detection of the extracted trend components was implemented using the Seasonal Hybrid Extreme Studentized deviate (S-H-ESD) method, which was built upon the Generalized ESD test for detecting anomalies. Machine learning was completed using the "AnomalyDetection" packages in R developed by Twitter (https://github.com/twitter/AnomalyDetection) [35]. With the packages, anomaly detection was achieved by employing time series decomposition and using robust statistical metrics, viz., median with ESD as well as piecewise approximation for high frequency time series data. The R package can detect both global and local anomaly values by setting the maximum anomaly probability among the all results. In this study, we assumed that the maximum anomaly month number may be induced during the implementation of water allocation. Thus, considering both positive and negative anomaly values, the maximum anomaly probability was set to 0.12 conservatively. The significant level was 95%.

ANOVA analysis
To evaluate the effects of water allocations, three groups were set: before water allocation, during water allocation, and after water allocation. Due to the different lengths of duration of water transformation, we selected monthly time series data groups to represent the before (Ahead group), during (Current group), and after water allocation groups (After group). Significant variations of NDVI categories under different groups at each buffer distance level at the three sites were systematically analyzed by one-way ANOVA, and multiple comparisons were also conducted among groups under the t-test. The analyses were performed in R (http://cran.R-project.org/).

PCA analysis
PCA (Principal Component Analysis) was performed to investigate relatively important variables of all NDVI categories to identify the major NDVI variables responding to water allocation. The "pca3d" package in R was used for the PCA. In addition, the suitability of PCA for our data was tested and confirmed by performing the matrix of cumulative variance contribution. In the research, only four key variables for both PC1 axis and PC2 axis were obtained deliberately, and 95% individuals  Table 2 shows descriptive statistics of NDVI categories for each gird at three buffer distance levels at three sampling sites during the recorded periods. The results showed similar grid numbers at each buffer distance level in all three sites. During the period without water allocation, the NDVI values, except NDVI min , were clearly larger than those during the water allocation period. Some negatives NDVI min values were observed during the period of without water allocation, which implies the threat of water scarcity to vegetation. In addition, the relatively smaller NDVI values along the water flow from the Angzh reservoir were compared with that of the other two sites, showing the natural spatial differences of water resources in the watershed, as well as the assisted proof of negative NDVI values at the Angzh site.

NDVI Trends and Anomaly Detection
The Robust STR could tolerate outliers well and provide details of NDVI trend components (marked as blue circles) at the Angzh (Fig.2), Wangk (Fig.3), and Xidayang (Fig.4) sites. A common was the stronger fluctuation during the last years of the recorded duration, especially towards the end of 2009. The different buffer distances resulted in different NDVI responses, with a huge fluctuation gradient from bf2 to bf6, especially in sites Angzh (Fig.2) and Xidayang (Fig.4). The marginal box plots show the trend data distribution and indicate differences of both spatial effects and distance effects. The box plots of all sites show relatively centralized data distribution for NDVI mean , NDVI max and NDVI most , and relatively higher values with increasing buffer distance, especially in NDVI max . In addition, a relatively narrower range of NDVI max and NDVI most was observed at the three buffer distance levels for all three sites, compared with the other two NDVI categories. The effects of water allocation are shown with red segments in the figures; no obvious NDVI trend was observed in response to water allocation. The simple linear regression of NDVI trend components indicated extremely significant decline trends (p<0.001) for all NDVI categories at all buffer distance levels in all sites, except NDVI min on bf2 and bf4 at Angzh (Fig.2c). The significant decline trends of all NDVI categories show the insufficient water resource in the watershed.
The anomaly detection for the trend components from STR showed consistent variations in NDVI max at the end of 2009 for all three buffer distance levels (Fig.5). In the anomaly detection, we assumed the maximum anomaly probability as 0.12, which means that the anomalies are 12% at most. The anomaly results showed global anomaly values of NDVI max , whereas local anomaly values appeared only at bf2 around 2000 in site Angzh (Fig.5Aa). Some scattered anomaly values also appeared in other NDVI categories, such as in NDVI mean at bf2 in site Xidayang (Fig.5C-a), NDVI mean at bf4 and bf6 in site Angzh (Figs.5A-bc), NDVI most at bf4 in site Angzh (Fig.5A-b), and so on.

Vegetation Ecological Responses and Buffer Distance Effects
The ecological responses of vegetation to water allocation in the three groups at different water transformation time in each site were analyzed using one-way ANOVA, comparisons were made between groups. As shown in Figures 6-8, extremely significant differences were observed between the groups at all buffer distance levels in all sites, except for that between the Ahead group and After group. At the beginning of the recorded water allocations (before July 2001) in sites Angzh and Wangk, NDVI values were relatively higher when receiving the transformed water than that in the same period before the water allocation (Figs. 6A, 6B, 7A). Moreover, NDVI values increased with the extension of buffer distance. However, at the end periods of the recorded water allocations (starting from 2009) opposite changes were detected in the comparisons between the Ahead group and Current group in Angzh site (Fig.6E). Among other comparisons of the three groups, the results showed nosignificant differences: Ahead group and Current group from June to July 2001 in Wangk site (p>0.05) (Fig.7B), Ahead group and After group from March to April 2005 in Angzh site (p>0.05) (Fig.6C), and Ahead group and After group from February to May 2005 in Angzh site (p>0.05) (Fig.8). The remaining comparatives groups showed extremely significant differences (p<0.01). Moreover, with the increase of buffer distance, the p value of t-test became relatively higher, which indicates the effects of distance on vegetation acquiring water. The results imply the benefit of vegetation restoration with water resource availability after the implementation of water transformation. In addition, we should emphasize on the importance of water allocation time. In the research, the response of vegetation to water acquisition was weak in spring (Figs. 6C,7E) and especially in summer, which is the water transformation time (Figs 6E,7B).

Relative important variables for ecological effects by PCA
A total of 12 variables in each site (combinations of four NDVI categories with three buffer distance levels) showed the ecological responses to water resource availability. The PCA could identify the relative important variables for their ecological sensitivity to water resource availability. As shown in Figure 9, NDVI min at bf4 and NDVI max at both bf4 and bf6 in Angzh are the most influential variables (Fig.9A). In Wangk, NDVI min at bf2 appears to be relatively more sensitive to water resources (Fig.9B). The variables of NDVI max at all buffer distance levels and NDVI min at bf4 display relatively higher sensitivity to water resources in Xidayang site (Fig.9C). The results reflect the effects of water transformation strategies aimed at vegetation restoration. When water allocation was implemented via the Wangk channel, riparian vegetation close to rivers was sensitive. While water transformation was accomplished by means of the Angzh channel, nonriparian vegetation far from rivers was sensitive to water resource.

Discussion
The analysis of descriptive statistics showed that all NDVI categories (except NDVI min ) had relatively high values in all sites during the period without water allocation. Omute et al. [23] revealed that NDVI could be a robust predictor for water level variations during low precipitation years around Lake Victoria. This result indicates the strong correlation of NDVI with water availability when natural water resources are not sufficient. Wang et al [36] revealed peak values of NDVI during the growing season when NDVI was highly correlated with precipitation on a monthly scale. Our results were similar only for NDVI min . In fact, water allocation is mostply implemented during the non-growing season when vegetation has less greenness representing the corresponding lower NDVI values. The negative NDVI value in the period without water allocation implies the threat of water scarcity to vegetation because of the relatively moistened conditions during the water allocation period. In addition, the negative NDVI value may represent more bare land rather than water bodies.
The NDVI trend by Robust STR detection revealed stronger fluctuation during the last years of the recorded duration, especially in 2009. Similar results were also reported by Wang et al [36]; RDA (redundancy analysis) results showed that hydrological conditions, especially water availability, had important effects on NDVI dynamics. As the economic and agriculture demand for water increased, water crisis occurred in all watersheds in the last years of the recorded data. Supporting evidence showed that relative small discharge flowed into the lake Baiyangdian in 2009, and evidence from the anomaly detection showed that global consistent variations of NDVI max occurred at the end of 2009. Our results also showed the effects of distance to water availability on vegetation with an obvious decline gradient from bf2 to bf6. Ling et al [28] reported significantly positive correlation of ecological water requirement with distance from rivers on both banks. The ecological water requirement of riparian vegetation can be satisfied over a stretch between 4.9 km and 10.0 km from the river [28]. Our results were similar with this result at least within the range of 6 km from rivers. The ANOVA results also indicated the effects of distance on vegetation response to water availability. Moreover, the p value of the t-test became relatively higher as the distance from the river was decreased, indicating the effects of distance on water acquisition by vegetation. Vegetation close to rivers could consume more water than that stored in soil during dry conditions compared with vegetation farther away from rivers.
Regional water resources have undergone substantial spatial and temporal variations since the implementation of water allocation. The intermittent frequency with different flow paths caused the spatial and temporal differences of vegetation distribution. Our results provide evidence for the spatial and temporal differences of vegetation responses to water availability. For example, NDVI values were relatively higher during the allocation period in Angzh beginning from 2000, while opposite changes occurred after 2008. Similar outcomes were also observed at different buffer distances. The above conclusions are also supported by the relative importance detection using PCA. For example, NDVI min at bf4 and NDVI max at both bf4 and bf6 in Angzh were the most prominent variables (Fig.9A), whereas NDVI min at bf2 had relatively higher sensitivity to water resource in Wangk (Fig. 9B). The PCA results also implied the effects of water allocation on vegetation at different buffer distances, which should be considered when the exact water flow path must be determined from multiple pathways with emphasis on large-scale vegetation restoration.
In addition, it should be noted that this study was somewhat limited by the raw and approximate hydrological data and the temporal and spatial resolution of the remote sensing data. The former data are intermittent records at an approximately month scale. The later data are representative of 10-day periods and we acknowledge that the temporal resolution contributes to uncertainty in the estimated significant changes in vegetation. Nevertheless, the extractions of NDVI should be robust because we used an unbiased method. In terms of spatial resolution, each 1-km 2 spatial resolution generally consisted of several different types of vegetation. It is important to acknowledge that the extracted NDVI values are representative of vegetation communities rather than individual species.

Conclusions
We analyzed the ecological effects of intermittent water allocation emphasizing on both temporal and spatial responses of the corresponding vegetation using remote sensing data and GIS-based buffer technology. The STR trend extraction and its corresponding linear regression and anomaly detection methods were used to determine temporal variations of vegetation under the effects of water allocation. ANOVA and PCA methods were used to identify the spatial responses of vegetation in different water flow paths and buffer distances. Our conclusions are as follows: (1) NDVI, except NDVI min , displayed higher values during the period without water allocation. The results imply implied water scarcity in the watershed and suggest an urgency for ecological restoration.
(2) The extremely significant decline trends (p<0.001) of all NDVI categories (except NDVI min ) at all buffer distance levels in all three sites at bf2 and bf4 in Angzh showed stronger fluctuations after 2008, as well as the decline gradient with buffer distance from rivers. The anomaly detection provided similar evidence of stronger NDVI fluctuations after 2008.
(3) The ecological response of vegetation to water allocation was significant, with regional vegetation coverage (p<0.01) showing a decline gradient of statistical p values along wide buffer distances. Our results provide evidences for the spatial and temporal differences of vegetation responses to water availability due to the intermittent frequency water allocation implemented via different river channels.
This study can be extended to other areas of the world where water diversion is being implemented or have been implemented. The data analysis in this study could be easily and conveniently replicated considering the convenience and high availability of global remote sensing data. The effects of water allocation on vegetation dynamics and wetland ecosystem restoration could also be conducted in other similar regions. Therefore, the implications of this work are potentially global. The findings presented herein are valuable for deepening our understanding of the effects of water diversion on regional vegetation restoration and can be used to develop a practical strategy for effective implementation of water allocation.