Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River

Quaternary International xxx (2015) 1e9

Contents lists available at ScienceDirect

Quaternary International journal homepage: www.elsevier.com/locate/quaint

Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River Ming Han a, b, Chengyi Zhao a, *, Gary Feng a, 1, Fengzhi Shi a a b

State Key Laboratory of Oasis Desert and Ecology, Xinjiang Institute of Ecology and Geography, CAS, Urumqi 830011, China University of Chinese Academy of Sciences, Beijing 100049, China

a r t i c l e i n f o

a b s t r a c t

Article history: Available online xxx

The responses of eco-hydrological system to anthropogenic and natural disturbances have attracted much attention in recent years. The coupling and simulating feedback between hydrological and ecological components have been realized in several recently developed eco-hydrological models. However, little research has been carried out to study the estimation of threshold parameters in the ecohydrological models. The aim of this paper is to infer the groundwater threshold parameter that connects the ecological and hydrological processes inside eco-hydrological models. A simplified vegetation dynamic model was set up, and used to simulate the vegetation dynamic along the lower reach of the Tarim River. Bayesian inference approach was applied with Markov Chain Monte Carlo sampling method used to infer the probability distribution of the groundwater threshold parameter. The result showed that the simplified model has the capacity to model the vegetation dynamic in the study area. The mean values of inferred groundwater threshold parameters for Yinsu, Kardayi, Alagan, and Yiganbjima are 4.99, 5.59, 5.74, and 5.85 m respectively, which are comparable with prior research, and showed significant spatial variability. The 90% confidence interval of the groundwater threshold parameter is larger than 1.27 m. This threshold parameter inference approach was suggested to be carried out before complex ecohydrological modeling is performed. © 2015 Elsevier Ltd and INQUA. All rights reserved.

Keywords: Bayesian inference Groundwater threshold Markov Chain Monte Carlo Eco-hydrological Vegetation dynamic

1. Introduction Riparian and wetland ecosystems has important environmental value, as they provide various ecosystem service such as mitigating floods, supporting wildlife habitats, and the most importantly, they prevent desertification in arid regions (Zhao et al., 2005; Chui et al., 2011; Sun et al., 2011; Shi et al., 2012). However, the riparian and wetland are threatened, because of human exploitation and climate change (Davis and Ogden, 1994; Chui et al., 2011; Van der Valk, 2012). The understanding of the riparian and wetland systems and their response to anthropogenic and natural disturbances are of foremost importance for management and restoration of riparian ecosystems (Wright and Chambers, 2002; Ridolfi et al., 2006). Integration of hydrology and ecology is needed, which is considered in the interdisciplinary research area of eco-hydrology (Crawford, 2000; Zalewski, 2002; Hannah et al., 2004).

* Corresponding author. E-mail addresses: [email protected] (C. Zhao), [email protected] (G. Feng). 1 Current address: Genetics and Precision Agriculture Research Unit, USDA-ARS, P. O. BOX 5367, 810 Highway 12 East, Mississippi State, MS 39762, USA.

Eco-hydrological models provide a potentially useful tool in characterizing hydrologicaleecological interactions. The prior ecohydrological modeling research could be characterized into three classes: (1)Some studies have added the vegetation water stress function into a hydrological model (Porporato et al., 2001; Rodriguez-Iturbe et al., 2001; Laio et al., 2001a, 2001b; Porporato et al., 2002; Brolsma and Bierkens, 2007; Lowry et al., 2011). Based on the assumption that ecosystems strive to minimize the vegetation stress, this approach was widely used to predict the vegetation pattern under climate change and human activities. Franz et al. (2010) concluded that parameters of water stress function played an important role in determining species distribution patterns, and the largest challenge that remained for this kind of modeling framework is the availability of species-specific parameters and the understanding of how those parameters varied in space; (2)The second approach incorporates the vegetation dynamic model into the hydrology model (Caylor et al., 2005; Ursino, 2005; Baudena et al., 2007; Liu, 2011; Scarsoglio et al., 2012). The environmental capacity, which is a variable inside the vegetation dynamic model, is the linkage between ecology and hydrology systems. A water stress function is included to accounts for the

http://dx.doi.org/10.1016/j.quaint.2015.02.035 1040-6182/© 2015 Elsevier Ltd and INQUA. All rights reserved.

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

2

M. Han et al. / Quaternary International xxx (2015) 1e9

dependence of environmental capacity on environment factors such as: the groundwater depth or soil moisture. This approach has been used to study the stability, resonance, and coherence in groundwateredependent ecosystems. Ridolfi et al. (2006) showed that the parameters in the water stress function and the shape of function had strong impact on the stability of eco-hydrological system; (3) The third approach couples a physically based vegetation growth model into the hydrology models, and the water stress function that influenced the actual root uptake is the linkage between two models (Wattenbach et al., 2005). In all of the above ecohydrological modeling approaches, the water stress function is the linkage between the ecology and the hydrology processes. This threshold of the water stress function is the core of the ecohydrological modeling. This threshold distinguishes the water use strategies among vegetation, which have strong impacts on the modeling result of vegetation competition, vegetation distributions, and hydrology processes. Further, this threshold is important for water resource management. However, little research has been carried out to study the estimation of threshold parameters in the water stress function. This kind of threshold has strong spatial variability, is heavily influenced by the local long term environment background, and is seldom measured directly (Guisan and Thuiller, 2005). Here, we introduce a Bayesian inference approach for threshold parameter estimation, which is a well-established and coherent methodology for statistical inference. A simplified vegetation dynamic model has been built for this purpose. The Tarim River, with a length of 1321 km, is located in the arid zone of northwestern China. Due to the severe abuse of water resources in the upper and middle reaches, the lower reach of Tarim river has dried up gradually, the groundwater table has dropped sharply, and the riparian vegetation has been degraded severely since the 1970s (Chen et al., 2010). The central and local governments of China have invested RMB 10.7 billion to restore and reconstruct the riparian ecosystem in the lower reach of the Tarim River. One of the key actions is water divergence from the Daxihaizi reservoir to the lower reach of the Tarim River. Until 2007, the water delivery has been carried out for 7 times. In support of such planning, research has been carried out. Research showed that the riparian ecosystem in the lower reach of the Tarim River is a groundwater dependent system, and the water conveyance considerably elevated the groundwater table and restored the riparian forest (Zhao et al., 2004; Chen et al., 2008, 2010; Hao et al., 2009; Yu et al., 2012). So the groundwater depth below which the system would suffer water stress is one of the key considerations when deciding the water convergence strategies. The water transfers enhance groundwater and allow riparian vegetation to develop at a much faster rate than under natural conditions (Liu et al., 2012). These significant interactions between ecological and hydrological processes in the Tarim River basin serve as a good experiment on eco-hydrological system evolution, and provide the possibility to apply a Bayesian approach to infer the groundwater threshold in this region. The main purposes of this paper are as follows: (1) Simulate the vegetation dynamics from 2000 to 2006 in four sections of lower reach of the Tarim river. (2) Apply Bayesian approach to infer threshold parameters in water stress function for the lower reach of the Tarim River. 2. Materials and methods 2.1. Model description To reduce the model parameters and more efficiently estimate the threshold parameter, our simplified vegetation dynamic model

only contained two parts: (1) vegetation dynamic model, (2)the connection between vegetation dynamic and groundwater depth. Instead of modeling the groundwatereriver interaction system, the groundwater observation was directly used as input. The Verhulst logistic model has been widely used as a basic equation to describe the vegetation dynamic in various ecohydrological models (Ridolfi et al., 2006; Muneepeerakul et al., 2008). Primarily, two variables are inside this model: Intrinsic growth rate b which accounts for the temporal response of vegetation, and the environment carrying capacity Vcc, which is the maximum amount of vegetation index sustainable with the available resource. The environment carrying capacity Vcc is the key variable that connects the environment factors such as groundwater depth to the vegetation dynamic. The increase or decrease of the vegetation index was defined by the difference between environment capacity and vegetation index status at a certain time. The vegetation dynamic model was defined as:


 dV V ¼bV 1 dt Vcc ðht Þ

(1)

Where: V is the vegetation dynamic index. ht is the groundwater depth (m) at time step t. The vegetation dynamic index has relatively wide context, and several indexes have been chosen to represent vegetation dynamic such as: biomass, vegetation coverage, Normalized Difference Vegetation Index (NDVI) (NoyMeir, 1975; Liu et al., 2012; Jeevarathinam et al., 2013). In ecological research, the NDVI derived from the satellite image was widely used as a tool to reflect the ecosystem status. For example, it has been used to study the woody plant response to drought (Lloret et al., 2007), and used to evaluate the post fire vegetation recovery processes (White et al., 1996). Recently, work has showed that the NDVI had significant response to the water convergence and could be used to represent the riparian vegetation recovery process (Wu and Tang, 2010; Sun et al., 2011). Here, the NDVI from MODIS product was used. An annual time step was used in the vegetation dynamic model. Yearly maximum NDVI in each year was used as vegetation index (V in eq (1)) to represent the vegetation dynamic, and seasonal averaged groundwater depth was used as model input. The impact of groundwater table on vegetation dynamic is difficult to quantify. When the groundwater table is too shallow, vegetation would suffer water stress, due to insufficient aeration of the root zone. If the groundwater table is too deep, vegetation can suffer water stress as well, due to lack of water supply. In most ecohydrological models, these effects are taken into account by a quadratic dependence of the environment carrying capacity Vcc, on the groundwater table depth, ht (Ridolfi et al., 2006; Scarsoglio et al., 2012). Due to the quadratic function properties, there is only one groundwater depth, at which the environment carrying capacity Vcc would reach its maximum value A. Conceptually, the suitable groundwater depth should be a range. Instead of defining the carrying capacity, Vcc as a quadratic function of ht, it is defined as piecewise functions of ht, as below:

Vcc ðht Þ ¼

8 > > > > > > < > > > > > > :

A

ht ha

ht < ha

A

hb  ht  ha

A k  ðht  hb Þ þ 1

ht > hb

(2)

Where: k is a shape coefficient; when groundwater depth is between ha and hb, the environment carrying capacity Vcc reach its maximum value A; when groundwater depth is smaller than ha, the

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

M. Han et al. / Quaternary International xxx (2015) 1e9

environment carrying capacity is linearly decreasing; when groundwater depth is larger than hb, the environment carrying capacity is decreasing as a hyperbolic function (Fig. 2).

  1 m m b2   s b 2 L ðDjbÞ ¼   lnð2pÞ   ln s 2 2 2 m X  ðDt  Vðxt ; bÞÞ2

3

(6)

t¼1

2.2. Formal Bayesian inference 2.2.1. Bayesian inference approach A general Bayesian statistical framework could be used to infer the unknown groundwater depth threshold in the vegetation dynamic model. Let D be an m dimensional vector containing the observed value of yearly maximum NDVI at time step t, t ¼ 1 … m. And b is the structural unknown parameters set in the vegetation dynamic model (b ¼ b, k, ha, hb, A). Let the m dimensional forcing data vector x (seasonal averaged groundwater depth) to be processed through vegetation dynamic model V(.):

D ¼ Vðx; bÞ þ ε

(3)

Where: V(x, b) is the vegetation response function described in section 2.1(eq. (1)), ε is the residual errors that represent errors caused by model structural, forcing data errors, and observation errors. In this study, model structural and forcing data errors were neglected and only observation error was considered. We wanted to infer the unknown parameter b from observation D. In the Bayesian approach, instead of assuming that b is a fixed parameter, we think it as a random variable. Our objective was to find the posterior probability distribution for unknown parameter set b, which captures the probabilistic beliefs about the parameters b in the light of the observed data D. According to Bayesian analysis framework, the posterior parameter distribution defined as (Kuczera and Parent, 1998; Gelman et al., 2013):

pðbjDÞ ¼ Z

pðDjbÞ  pðbÞ pðDjbÞ  pðbÞdðbÞ

fpðDjbÞ  pðbÞ

(4)

Where: pðbjDÞ is the posterior parameter distribution, is the likelihood function, p(b) is the prior probability density distribution.

2.2.2. Likelihood function The likelihood is a statistical measure of the quality of the fit between model result and the observation. It describes the probability of obtaining the observations D given the model parameter b (Hartig et al., 2012; Gelman et al., 2013). Various likelihood functions have been used and tested for different purposes (Freer et al., 1996; Vrugt et al., 2003; Gelman et al., 2010). If we assume that the error residuals to be uncorrelated, Gaussian distributed with constant variance s2 and zero mean, the likelihood function can be defined as (Vrugt et al., 2003; Hartig et al., 2012; Vrugt and Sadegh, 2013):

pðDjbÞ ¼

  1 1 b 2  ðDt  Vðxt ; bÞÞ2 pffiffiffiffiffiffiffiffiffiffiffiffi  exp   s 2 b2 2p s t¼1 m Y

(5)

b is an estimate of the standard deviation of the meaWhere: s b was specified as 0.025, based surement error. The value of the s on knowledge of the MODIS NDVI product measurement errors (NASA, 2012; Vrugt and Sadegh, 2013). Dt is the observation of yearly maximum NDVI at year t. m is the number of the measurement. For reasons of calculation simplicity and numerical stability, it is often consider the log likelihood function, which defined as:

2.2.3. Sampling method In this case, it was difficult to calculate the posterior distribution pðbjDÞ directly, since we do not know the normalizing constant. The Metropolis-Hastings algorithm, which is the basic building block of classical Markov chain Monte Carlo (MCMC) methodology, was used to estimate pðbjDÞ. The key idea of the MCMC is building a Markov chain by generating the sequence of samples {b1…bN}. Due to its ergodic properties, it can be shown that it converges to the stationary distribution, pðbjDÞ. This method already has received considerable attention in the last decade in the Bayesian statistics literature (Kuczera and Parent, 1998; Gelman et al., 2013). The “random-walk” Metropolis algorithm, which was introduced by Tierney (1994), was used to explore the posterior distributions in this study. It is implemented as follows: Step 1: For the ith iteration, randomly select the increment Z from the increment distribution T ðÞ, which was a Multi-normal distribution on the parameter space, with mean zero, and covariance matrix I, where I is an identity matrix:

Z)T ðÞ

(7)

Then the new candidate bc in time step i was defined as:

bc ¼ bi1 þ scale  Z

(8)

Where: scale is the user specified scale factor, which is very sensitive to acceptance rate. The only criterion known for choosing sensible scaling is to adjust it so that the acceptance rate is about 20% (Geyer and Johnson, 2012). In this study a scale ¼ 0.05 was used. Step 2: Evaluate the acceptance probability (u) from the ratio:

pðbc jDÞ u ¼ min 1; pðbi1 jDÞ

(9)

Step 3: Randomly sample from a uniform distribution over the interval 0 to 1

R)Uð0; 1Þ

(10)

If R  u, then jump to bc, that is bi ¼ bc. However if R > u, then remain at current position, that is bi ¼ bi1. This aforementioned Metropolis-Hastings algorithm was included in the R package “mcmc” (Geyer and Johnson, 2012). After we generated samples {b1…bN} following the postal destitution of pðbjDÞ, the marginal probability distribution of groundwater threshold parameter hb, pðhb jDÞ was estimated by using relative frequency of hb in every sub-interval of hb. The parameter space of hb, ranging from 2 to 8 m, was divided into 80 sub-regions with 0.1 m interval. 2.3. Data collection and model setup 2.3.1. Study area The study area is located between the Dxihaizi reservoir and Taitema Lake in lower reach of the Tarim River. It is one of the most arid areas in China, with an average annual precipitation less than 50 mm, the potential annual evaporation ranging between 2500 mm and 3000 mm (Fig. 1) (Chen et al., 2010). Vegetation is

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

4

M. Han et al. / Quaternary International xxx (2015) 1e9

distributed along the river belt. Common flora consists of trees, shrubs, and herbs. Trees include Populus euphratica, shrubs include Tamarix spp and Lycium ruthenicum, and herbs include Phragmites communis (Chen et al., 2010). In order to evaluate the ecological restoration since 2000, nine observation sites were established for monitoring groundwater levels and vegetation dynamics (Fig. 1). The distance between the two neighboring sections of the first 6 sections is about 20 km, and the distance between the last three sections is about 40 km (Chen et al., 2010). All sections are at 90 to the channel. Three to eight groundwater wells were installed along the centerline of each transaction in distances of 50 m, 150 m, 350 m, 450 m, 700 m, and 1000 m respectively (Chen et al., 2010). The groundwater table was measured using the conductance method and we generally took measurements every 1 or 2 months, but every 10 days during the water conveyance (Chen et al., 2010). The vegetation dynamic model and the Bayesian inference framework was applied at monitoring sections Yinsu, Kaerdayi, Alagan, and Yiganbujima, respectively, which are spatially distributed from upper part of study area to lower part of study area, from 2001 to 2006. In each section, the modeling site is at a groundwater observation well that was 150 m away from the river bank. 2.3.2. Vegetation index (Yearly maximum NDVI) The NDVI data of MODIS has been widely used in ecological studies. In the current work, NDVI measurement from MODIS/Terra MOD13Q1 product from 2000 to 2006 was used. The spatial resolution of this product is 250 m and the temporal resolution is 16 days. For each section, the time series of MODIS cells that correspond to groundwater observations wells which is 150 m away from river bank was extracted, and then yearly maximum NDVI value from the time series was calculated and used to represent the vegetation dynamic. The yearly maximum NDVI of 2000 was used as initial

Fig. 2. Relationship between environment carrying capacity Vcc and groundwater depth.

condition and the remaining maximum NDVI time series were used to calibrate vegetation dynamic model and inferred threshold parameter hb.

2.3.3. Groundwater data In this groundwater-dependent ecosystem, the depth of the groundwater was a strong predictor of plant frequency. In each section, groundwater depth measurements from an observation well that is 150 m away from the river bank were used. The yearly average groundwater depths in vegetation growth period covering from July to September were calculated from 2001 to 2006. This yearly data set was used as model input for the vegetation dynamic model.

2.3.4. Prior information and convergence diagnostic Totally, there were five parameters (b, k, ha, hb, A) inside the above simplified vegetation dynamic model. However, during the simulation, the threshold ha had no effect during the inference process, mainly because the groundwater depth was elevated from about 9 m to about 4 m in the study area from 2000 to 2006, and as a result, the case of h < ha (see eq. (2)) did not occur. Thus, only four parameters (b, k, hb, A) were inferred in this study. A uniform prior [p(b) ¼ 1] over the range of parameters was used, and prior ranges of the parameters were defined in Table 1. The most critical issue of the implementing MCMC method is to determine whether convergence has been achieved. To assess the convergence of Markov Chain, the R statistic (Gelman and Rubin, 1992) was used as a quantitative convergence diagnostic, which they called the scale reduction score. The R statistic is a measure of the between-chain and within-chain variances. If the multiple chains are not mixing properly in the parameter space the R statistic will be relatively high, and as the R statistic close to 1 for each of the parameter indicating convergence. However, because a score of unity is difficult to achieve, Gelman and Rubin (1992) recommended using values less than 1.2 to declare convergence to a stationary distribution.

Table 1 Prior Uncertainty ranges of the eco-hydrology model parameters.

Fig. 1. Study area and monitoring sections along lower reaches of Tarim River.

Max. Min.

A [e]

hb [m]

b [e]

k [e]

0.2 0.1

7 2

0.8 0.01

3.000 0.010

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

M. Han et al. / Quaternary International xxx (2015) 1e9

In all Metropolis-Hastings applications, the following procedure was used: 5 parallel chains with different initial starting pointing which were over dispersed in the parameter space were set up. Each chain had 2,000,000 samples. After confirming convergence, the final joint posterior distribution of pðbjDÞ, was obtained by combining all 5 chains samples after discarding the initial burn-in period (the first half of the Markov chain). In total, 5,000,000 samples were used for estimating the joint posterior distribution of pðbjDÞ in each section. 3. Results 3.1. Convergence of MCMC To evaluate the convergence of the Markov chains, the R statistic value for parameters was calculated in each section. All R statistic values were smaller than 1.2, indicating Markov chains were converged to the stationary posterior distribution. The R statistic values of parameter hb were 1.08, 1.04, 1.03, and 1.11 in Yinsu, Karedayi, Alagan, and Yinganbujima respectively. The evolution of the R statistic value and the trace of Markov chains samples for hb were further plotted in Kaerdayi as an example to illustrate the convergence. For Kaerdayi, the evolution of the R statistic values for all parameters are illustrated in Fig. 3A. The R statistic values were larger than 30 at the beginning, because of random initializing of the starting points of sequences. It decreased quickly, and after 10,000 simulations all R statistic values were smaller than 3. This part was not mapped to more clearly trace the R statistic value. Beginning with the 10,000th simulation, the R statistic value continued increasing and decreasing intermittently. Finally, after approximately 300,000 simulations, the R statistic values of parallel sequences in the plot were smaller than 1.2, suggesting that the all parameters had converged to a stationary posterior distribution.

5

Three of the five parallel sequences of parameter hb are illustrated in Fig. 3B. The plot shows that at early stages during the evolution, the individual sequences tend to occupy different regions of the posterior surface before 10,000 simulations. This low mixing of the paths was associated with a relatively high value for the R statistic value, indicating poor convergence (Fig. 3A). At a later stage, after 200,000 simulations, all of the individual sequences have been able to fully explore the posterior target distribution and showed a high mixing path, thereby resulting in a R statistic value smaller than 1.2, indicating convergence to a stationary distribution. The same analysis has been conducted for other sections. All R statistic values were smaller than 1.2 after 1,000,000 simulations. The samples taken from 1,000,000 simulations to 2,000,000 simulations could provide a good approximation to the posterior distribution, and are sufficiently accurate to support the analyses presented in this paper (Gelman and Rubin, 1992). 3.2. Parameter estimation Fig. 4 presents the posterior marginal probability density distribution for hb inferred for each section using the samples generated with the Metropolis-Hastings algorithm. The sharp and peaked distribution of parameter hb at Kaerdayi and Yinsu is well identifiable, while flat distributions for Alagan and Yiganbujima indicated more parameter uncertainty. This was confirmed by the standard deviation of parameter samples given in Table 2, which was the summary statistic for parameter hb at all sections. From Fig. 4 and Table 2, the mean values of inferred threshold parameters hb for Yinsu, Kardayi, Alagan, and Yiganbujima were 4.99, 5.59, 5.74, and 5.85 m respectively. Parameter hb showed significant spatial variations. The maximum difference of mean value of hb among sections was 0.86 m (between Yinsu and Yiganbujima). The possibility that Yinsu and section Yiganbjima had the same value of parameter hb was 0.025. From the upper reach (Yinsu) to the lower

Fig. 3. Evolution of the Gelman and Rubin scale-reduction factor for the parameters in section Kaerdayi (A), Markov Chain Monte Carlo samples of hb generated in three parallel sequences using Metropolis-Hastings algorithm (B).

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

6

M. Han et al. / Quaternary International xxx (2015) 1e9

Fig. 4. Marginal posterior probability distributions of the parameter hb for each section.

reach (Yinganbjima), the mean value of hb showed an increasing trend. This may be caused by long term spatial distribution of groundwater depth among sections during the last 50 years, which was also increasing from upper reach (Yinsu) to lower reach (Alagan) (Chen et al., 2006). This reveals the vegetation adaptation strategies depend on the environment factor: the deeper the long term groundwater depth, the deeper the root system the vegetation would developed, and the larger the threshold parameter hb. A similar result has been reported in the control experiment for P. euphratica and Tamarix ramosissima seedlings (Li et al., 2013). Li et al. (2013) concluded that Root/shoot ratios and maximum root depth of both species tended to increase with increasing groundwater depth.

3.3. Model prediction uncertainty Finally, in order to verify the model capability to simulate the vegetation dynamic in these sections, the result of the Bayesian inference was translated into estimates of the model prediction

Table 2 Marginal posterior Mean (Mean), Standard Deviation (SD), Coefficient of Variation (CV), and 90% confidence interval for threshold parameter hb (m) in each section. Sections

Mean

SD

CV

P ¼ 0.05

P ¼ 0.95

Yinsu Kardayi Alagan Yiganbjima

4.99 5.59 5.74 5.85

0.414 0.305 0.581 0.644

0.083 0.055 0.101 0.110

4.24 5.08 4.69 4.64

5.59 6.06 6.64 6.83

uncertainty. For each section, 100,000 parameter sets were obtained from 5,000,000 samples by taking every 50. The model was run and the value of likelihood between simulation and observation was calculated for each of 100,000 thinned parameter sets. The 90% confidence intervals of model uncertainty were estimated by sorting the corresponding likelihood values (Jin et al., 2010). The observed NDVI and the 90% model prediction uncertainty interval are shown in Fig. 5. For all sections, the model prediction uncertainty interval included most of the observations, and captured the NDVI dynamic trends, suggesting that the model structure was suitable for simulating vegetation dynamic in our case study. The 90% confidence interval for parameter hb was large (Table 2), and the minimum 90% confidence interval was 1.27 m in Yinsu. Parameter hb inside the 90% confidence interval could give an equivalent simulation result (Fig. 5). We suggest that before complex eco-hydrological real system modeling, this approach should be applied to get insight concerning hb. 4. Discussion Most of the eco-hydrological modeling is “data hungry”, as numerous hard-to- measure physical parameters are laborious to obtain. Furthermore, model parameters often have strong spatial variability (Guisan and Thuiller, 2005). But on the other hand, many existing data, such as species coverage, remote sensing product (NDVI) and species distribution database, cannot be used because they describe variables that are modeled on the basis of multiple, interacting processes and therefore depend on the physical parameters (Hartig et al., 2012). Inverse modeling is a bridge between

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

M. Han et al. / Quaternary International xxx (2015) 1e9

7

Fig. 5. The simplified eco-hydrology model prediction uncertainty derived from MCMC result. The lighter shade region denotes the prediction 90% confidence interval, while the black dots corresponding to observed NDVI data.

observation and unknown parameters. Inverse modeling approaches could be classified into two classes: direct inverse modeling and statistical. Normal direct inverse modeling defines an objective function to measure the goodness of fit, and finds the best parameter set via optimization. Compared with this approach, the Bayesian inference statistical approach considers the unknown parameter as a random variable as well and captures the probabilistic beliefs about the unknown parameters in the light of the observed data. Due to measurement errors, model simulations never perfectly represent the system or exactly fit the data. There may be many parameter sets within a model structure that are equally acceptable as simulators of the real system, and these may often come from very different region in the parameter space (Freer et al., 1996). As Kuczera and Parent (1998) noted: “no serious hydrologist should be naive enough to rely on a uniquely determined value for model parameters”. The posterior distribution of parameters provides more than a simple point estimate or even a confidence interval. Instead, it describes the relative probabilities of various portions of the parameter space and provides information on the full correlation structure, which can be difficult to access in other analysis. Compared with other statistical inversing modeling approach, Bayes' theorem offers the possibility of merging some

independent knowledge about parameters into the inference process. The posterior distribution of parameters generated from Bayesian theory is easy to update when another source of information is available. Due to the above advantage, Bayesian inference has been widely applied in hydrology and ecology studies (ter Braak and Etienne, 2003; Rivot et al., 2004; Ershadi et al., 2013). Most previous research inferred the maximum suitable groundwater depth though studying the relationship between groundwater observation and vegetation diversity, using the groundwater data and field investigation data at all transects (Hao et al., 2009, 2010; Ma et al., 2011). According to their result the suitable groundwater depth was between 4 and 6 m. Our result shown in Table 2 was comparable with their research. However, the meaning of the threshold estimated by Hao et al. (2010) is the groundwater depth under which the vegetation has the maximum possibility of occurrence, while the meaning of our result is the groundwater depth below which the vegetation would suffer water stress. The advantage of the method in this paper is that: 1) Beside estimation of threshold parameters, our method could model the vegetation recovery process during 2000e2006 and the impact of groundwater depth on the vegetation dynamic. At the same time, our approach could also help the authorities to make a water

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

8

M. Han et al. / Quaternary International xxx (2015) 1e9

convergence plan. 2) Due to the complexity of the natural system, it is impossible for us to get the true value of the groundwater threshold. What we can do is to find ways to make the estimated threshold value closer to the true value. The result of this paper could be updated and refined, when other information such as vegetation coverage or longer time series data are available. 3) This approach could reveal the spatial variation of parameters. The inference approach applied in this paper could potentially combine spatial distributed groundwater model simulation results with spatial distributed observations such as remote sensing products, which would provide more insight of the spatial variation of the threshold parameters. Concerning the water management, our result of threshold parameter hb could be used for vegetation restoration in the lower reach of the Tarim River. However,t the suitable hydrologic environment that needed for riparian forest succession was not studied here. P. euphratica, which is the dominant vegetation along the riparian forest, shows sexual reproduction and asexual reproduction. The sexual reproduction relies on flooding events, and seed release is often synchronized with timing of flood events, which could provide ideal germination and initial growth conditions. We would suggest that a flood water influx should be applied in certain years. The amount of water required to create flooding events and the corresponding possible inundation region needs to be further studied. Wiehle et al. (2009) concluded that the root suckering of P. euphratica is a highly successful strategy to overcome the limited spatial and temporal gaps for seed germination and extend stands along the river. However, the water condition that may improve P. euphratica root suckers is still unclear, which is another key issue for sustainable riparian forest restoration. 5. Conclusions The threshold parameter, which regulates the vegetation responses to environments change, has drawn less attention in previous hydrology modeling, due to the complexity of the hydrology system and small impact of the threshold parameter on hydrology modeling results. With the development of the eco-hydrological modeling, the threshold parameter has shown its importance recently. It has significant impact on species patterns prediction, stability of eco-hydrological systems, and competition between species (Ridolfi et al., 2006; Franz et al., 2010). However, little work has been carried out to study the estimation and the uncertainty of the threshold parameter in eco-hydrological modeling. We applied a Bayesian inference approach to estimate the groundwater threshold for arid riparian vegetation in the lower reach of the Tarim River with a simplified vegetation dynamic model. The result showed that the simplified model had the capacity to model the vegetation dynamic in the study area. The mean values of inferred groundwater threshold hb for Yinsu, Kardayi, Alagan, and Yiganbjima were 4.99, 5.59, 5.74, and 5.85 m respectively, which were comparable with prior research, and showed significant spatial variability. The 95% confidence interval for parameter hb was large, and the minimum 90% confidence interval was 1.27 m in Yinsu. This parameter inference approach should be carried out before complex eco-hydrological modeling. Acknowledgements The study was supported by the National Natural Science Foundation (4117095, U14003281), National Project (2013BAC10B01) and National 973 Project (2013CB429905) of China. The authors are grateful to Dr Pengnian Yang of Xinjiang Agriculture University, China and Keilholz Patrick of Technical

University Munich, Germany, for their help in obtaining the groundwater observations. References Baudena, M., Boni, G., Ferraris, L., Von Hardenberg, J., Provenzale, A., 2007. Vegetation response to rainfall intermittency in drylands: results from a simple ecohydrological box model. Advances in Water Resources 30, 1320e1328. Brolsma, R.J., Bierkens, M.F.P., 2007. Groundwateresoil waterevegetation dynamics in a temperate forest ecosystem along a slope. Water Resources Research 43, W01414. Caylor, K.K., Manfreda, S., Rodriguez-Iturbe, I., 2005. On the coupled geomorphological and ecohydrological organization of river basins. Advances in Water Resources 28, 69e86. Chen, Y.N., Chen, Y.P., Xu, C.C., Ye, Z.X., Li, Z.Q., Zhu, C.G., Ma, X.D., 2010. Effects of ecological water conveyance on groundwater dynamics and riparian vegetation in the lower reaches of Tarim River, China. Hydrological Processes 24, 170e177. Chen, Y.N., Pang, Z.H., Chen, Y.P., Li, W.H., Xu, C.C., Hao, X.M., Huang, X., Huang, T.M., Ye, Z.X., 2008. Response of riparian vegetation to water-table changes in the lower reaches of Tarim River, Xinjiang Uygur, China. Hydrogeological Journal 16, 1371e1379. Chen, Y.N., Zilliacus, H., Li, W.H., Zhang, H.F., Chen, Y.P., 2006. Ground-water level affects plant species diversity along the lower reaches of the Tarim River, Western China. Journal of Arid Environments 66, 231e246. Chui, T.F.M., Low, S.Y., Liong, S.Y., 2011. An ecohydrological model for studying groundwaterevegetation interactions in wetlands. Journal of Hydrology 409, 291e304. Crawford, R., 2000. Eco-hydrology: plants and water in terrestrial and aquatic environments. Journal of Ecology 88, 1095e1096. Davis, S., Ogden, J.C., 1994. Everglades: the Ecosystem and Its Restoration. CRC Press. Ershadi, A., McCabe, M.F., Evans, J.P., Mariethoz, G., Kavetski, D., 2013. A Bayesian analysis of sensible heat flux estimation: quantifying uncertainty in meteorological forcing to improve model prediction. Water Resources Research 49, 2343e2358. Franz, T.E., Caylor, K.K., Nordbotten, J.M., Rodríguez-Iturbe, I., Celia, M.A., 2010. An ecohydrological approach to predicting regional woody species distribution patterns in dryland ecosystems. Advances in Water Resources 33, 215e230. Freer, J., Beven, K., Ambroise, B., 1996. Bayesian estimation of uncertainty in runoff prediction and the value of data: an application of the GLUE approach. Water Resources Research 32, 2161e2173. Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B., 2013. Bayesian Data Analysis. CRC Press. Gelman, A., Leenen, I., Mechelen, I.V., Boeck, P.D., Poblome, J., 2010. Bridges between deterministic and probabilistic models for binary data. Journal of the Royal Statistical Society B 7, 187e209. Gelman, A., Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences. Statistics 457e472. Geyer, C., Johnson, L.T., 2012. mcmc: Markov Chain Monte Carlo. http://cran.rproject.org/web/packages/mcmc/index.html. Guisan, A., Thuiller, W., 2005. Predicting species distribution: offering more than simple habitat models. Ecological Letters 8, 993e1009. Hannah, D.M., Wood, P.J., Sadler, J.P., 2004. Ecohydrology and hydroecology: a ‘new paradigm’? Hydrological Processes 18, 3439e3445. Hao, X.M., Chen, Y.N., Li, W.H., 2009. Indicating appropriate groundwater tables for desert river-bank forest at the Tarim River, Xinjiang, China. Environmental Monitoring and Assessment 152, 167e177. Hao, X.M., Li, W.H., Huang, X., Zhu, C.G., Ma, J.X., 2010. Assessment of the groundwater threshold of desert riparian forest vegetation along the middle and lower reaches of the Tarim River, China. Hydrological Processes 24, 178e186. Hartig, F., Dyke, J., Hickler, T., Higgins, S.I., O'Hara, R.B., Scheiter, S., Huth, A., 2012. Connecting dynamic vegetation models to dataean inverse perspective. Journal of Biogeography 39, 2240e2252. Jeevarathinam, C., Rajasekar, S., Sanjuan, M.A., 2013. Vibrational resonance in groundwater-dependent plant ecosystems. Ecological Complexity 15, 33e42. Jin, X., Xu, C.-Y., Zhang, Q., Singh, V.P., 2010. Parameter and modeling uncertainty simulated by GLUE and a formal Bayesian method for a conceptual hydrological model. Journal of Hydrology 383, 147e155. Kuczera, G., Parent, E., 1998. Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. Journal of Hydrology 211, 69e85. Laio, F., Porporato, A., Fernandez-Illescas, C., Rodriguez-Iturbe, I., 2001a. Plants in water-controlled ecosystems: active role in hydrologic processes and response to water stress: IV. Discussion of real cases. Advances in Water Resources 24, 745e762. Laio, F., Porporato, A., Ridolfi, L., Rodriguez-Iturbe, I., 2001b. Plants in watercontrolled ecosystems: active role in hydrologic processes and response to water stress: II. Probabilistic soil moisture dynamics. Advances in Water Resources 24, 707e723. Li, J., Yu, B., Zhao, C., Nowak, R.S., Zhao, Z., Sheng, Y., Li, J., 2013. Physiological and morphological responses of Tamarix ramosissima and Populus euphratica to altered groundwater availability. Tree Physiology 33, 57e68.

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035

M. Han et al. / Quaternary International xxx (2015) 1e9 Liu, D., Tian, F., Hu, H., Lin, M., Cong, Z., 2012. Ecohydrological evolution model on riparian vegetation in hyperarid regions and its validation in the lower reach of Tarim River. Hydrological Processes 26, 2049e2060. Liu, H.-H., 2011. Impact of climate change on groundwater recharge in dry areas: an ecohydrology approach. Journal of Hydrology 407, 175e183. Lloret, F., Lobo, A., Estevan, H., Maisongrande, P., Vayreda, J., Terradas, J., 2007. Woody plant richness and NDVI response to drought events in Catalonian (northeastern Spain) forests. Ecology 88, 2270e2279. Lowry, C.S., Loheide, S.P., Moore, C.E., Lundquist, J.D., 2011. Groundwater controls on vegetation composition and patterning in mountain meadows. Water Resources Research 47. Ma, X.D., Chen, Y.N., Zhu, C.G., Li, W.H., 2011. The variation in soil moisture and the appropriate groundwater table for desert riparian forest along the Lower Tarim River. Journal of Geographical Science 21, 150e162. Muneepeerakul, C.P., Miralles-Wilhelm, F., Tamea, S., Rinaldo, A., RodriguezIturbe, I., 2008. Coupled hydrologic and vegetation dynamics in wetland ecosystems. Water Resources Research 44. NASA, 2012. http://landval.gsfc.nasa.gov/. Noy-Meir, I., 1975. Stability of grazing systems: an application of predator-prey graphs. Journal of Ecology 459e481. Porporato, A., D'odorico, P., Laio, F., Ridolfi, L., Rodriguez-Iturbe, I., 2002. Ecohydrology of water-controlled ecosystems. Advances in Water Resources 25, 1335e1348. Porporato, A., Laio, F., Ridolfi, L., Rodriguez-Iturbe, I., 2001. Plants in water-controlled ecosystems: active role in hydrologic processes and response to water stress: III. Vegetation water stress. Advances in Water Resources 24, 725e744. Ridolfi, L., D'Odorico, P., Laio, F., 2006. Effect of vegetationewater table feedbacks on the stability and resilience of plant ecosystems. Water Resources Research 42. Rivot, E., Prevost, E., Parent, E., Bagliniere, J.L., 2004. A Bayesian state-space modelling framework for fitting a salmon stage-structured population dynamic model to multiple time series of field data. Ecological Modelling 179, 463e485. Rodriguez-Iturbe, I., Porporato, A., Laio, F., Ridolfi, L., 2001. Plants in water-controlled ecosystems: active role in hydrologic processes and response to water stress: I. Scope and general outline. Advances in Water Resources 24, 695e705. Scarsoglio, S., D'Odorico, P., Laio, F., Ridolfi, L., 2012. Spatio-temporal stochastic resonance induces patterns in wetland vegetation dynamics. Ecological Complexity 10, 93e101. Sun, Z., Chang, N.B., Opp, C., Hennig, T., 2011. Evaluation of ecological restoration through vegetation patterns in the lower Tarim River, China with MODIS NDVI data. Ecological Informatics 6, 156e163. Sun, D., Zhao, C., Wei, H., Peng, D., 2011. Simulation of the relationship between land use and groundwater level in Tailan River basin, Xinjiang, China. Quaternary International 244, 254e263.

9

Shi, F., Zhao, C., Sun, D., Peng, D., Han, M., 2012. Conjunctive use of surface and groundwater in central Asia area: a case study of the Tailan River Basin. Stoch Stochastic Environmental Research and Risk Assessment 26, 961e970. ter Braak, C.J.F., Etienne, R.S., 2003. Improved Bayesian analysis of metapopulation data with an application to a tree frog metapopulation. Ecology 84, 231e241. Tierney, L., 1994. Markov Chains for exploring posterior distributions. Annals of Statistics 22, 1701e1762. Ursino, N., 2005. The influence of soil properties on the formation of unstable vegetation patterns on hillsides of semiarid catchments. Advances in Water Resources 28, 956e963. Van der Valk, A.G., 2012. The Biology of Freshwater Wetlands. Oxford University Press. Vrugt, J.A., Gupta, H.V., Bouten, W., Sorooshian, S., 2003. A Shuffled Complex Evolution Metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters. Water Resources Research 39. Vrugt, J.A., Sadegh, M., 2013. Toward diagnostic model calibration and evaluation: approximate Bayesian computation. Water Resources Research 49, 4335e4345. Wattenbach, M., Hattermann, F., Weng, R., Wechsung, F., Krysanova, V., Badeck, F., 2005. A simplified approach to implement forest eco-hydrological properties in regional hydrological modelling. Ecological Modelling 187, 40e59. White, J.D., Ryan, K.C., Key, C.C., Running, S.W., 1996. Remote sensing of forest fire severity and vegetation recovery. International Journal of Wildland Fire 6, 125e136. Wiehle, M., Eusemann, P., Thevs, N., Schnittler, M., 2009. Root suckering patterns in Populus euphratica (Euphrates poplar, Salicaceae). Trees 23, 991e1001. Wright, J.M., Chambers, J.C., 2002. Restoring riparian meadows currently dominated by Artemisa using alternative state concepts-above-ground vegetation response. Applied Vegetation Science 5, 237e246. Wu, J., Tang, D., 2010. The influence of water conveyances on restoration of vegetation to the lower reaches of Tarim River. Environmental Earth Sciences 59, 967e975. Yu, P.J., Xu, H.L., Ye, M., Liu, S.W., Gong, J.J., An, H.Y., Fu, J.Y., 2012. Effects of ecological water conveyance on the ring increments of Populus euphratica in the lower reaches of Tarim River. Journal of Forest Research 17, 413e420. Zalewski, M., 2002. Ecohydrologydthe use of ecological and hydrological processes for sustainable management of water resources (EcohydrologiedLa prise en cologiques et hydrologiques pour la gestion durable des compte de processus e ressources en eau). Hydrological Sciences Journal 47, 823e832. Zhao, C., Wang, Y., Chen, X., Li, B., 2005. Simulation of the effects of groundwater level on vegetation change by combining FEFLOW software. Ecological Modelling 187, 341e351. Zhao, C., Wang, Y., Hu, S., Li, Y., 2004. Effects of spatial variability on estimation of evapotranspiration in the continental river basin. Journal of Arid Environments 56, 373e382.

Please cite this article in press as: Han, M., et al., Bayesian inference of the groundwater depth threshold in a vegetation dynamic model: A case study, lower reach, Tarim River, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.02.035